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13.  ABSTRACT  (Mwmwn 200 wordf) 

The  aim  of  this  thesis  is  to  calculate  and  discuss  some  of  the  properties  of 
probability  density  functions  of  bearing  angle  conditioned  on  the  data  received  by 
an  array  of  sensors. 

The  development  of  the  thesis  goes  through  three  major  stages.  The  first  stage 
of  the  development  is  the  theoretical  derivation  of  the  probability  density 
function.  The  second  stage  of  the  development  concerns  the  calculation  of  the 
density  function  with  examples.  This  development  includes  a  detailed  discussion 
of  the  simulation  used  to  produce  the  data  on  which  the  density  functions  are 
conditioned  and  the  code  written  to  do  the  actual  computations.  The  third  stage 
of  the  development  is  the  analysis  of  estimates  which  are  made  using  the 
calculated  density  functions.  This  stage  includes  comparison  of  the  results 
obtained  using  the  calculated  density  functions  and  a  better  known  performance 
measure:  the  Cramer-Rao  bound. 
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Abstract 


The  aim  of  this  thesis  is  to  calculate  and  discuss  some  of  the  properties  of  probability 
density  functions  of  bearing  angle  conditioned  on  the  data  received  by  an  array  of  sensors. 


The  development  of  the  thesis  goes  through  three  major  stages.  The  first  stage  of  the 
development  is  the  theoretical  derivation  of  the  probability  density  function.  The  second 
stage  of  the  development  concerns  the  calculation  of  the  density  function  with  examples. 
This  development  includes  a  detailed  discussion  of  the  simulation  used  to  produce  the  data 
on  which  the  density  functions  are  conditioned  and  the  code  written  to  do  the  actual 
computations.  The  third  stage  of  the  development  is  the  analysis  of  estimates  which  are 
made  using  the  calculated  density  functions.  This  stage  includes  comparison  of  the  results 
obtained  using  the  calculated  density  functions  and  a  better  known  performance  measure: 
the  Cramer-Rao  bound. 

The  major  results  obtained  are  as  follows.  The  first  major  result  is  that  the  density 
functions  of  bearing  conditioned  on  the  actual  data  received  at  an  array  of  sensors  are 
skewed  as  bearing  angle  increases  from  the  0°  array  centerline.  Better  estimates  are 
obtained  for  more  samples,  better  sampling  resolution,  or  higher  SNR.  The  final  result 
is  that  the  error  variances  of  estimates  made  using  the  density  functions  calculated 
decrease  with  increasing  sample  number,  SNR,  and  sampling  resolution. 
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Chapter  1  Outline  and  Goals 


1.1  Introduction 

The  work  in  this  thesis  is  developed  in  order  to  calculate  and  explore  the  properties  of 
probability  density  functions  of  the  bearing  angle  to  a  source  conditioned  on  data 
sampled  at  the  outputs  of  an  array  of  sensors.  The  thesis  begins  with  the  theoretical 
background  for  the  computation  of  the  density  functions.  This  includes  a  close 
examination  of  the  covariance  matrix  of  the  data.  The  development  then  turns  to  the 
actual  computation  of  the  density  functions,  the  code  used  to  construct  the  covariance 
matrix  and  the  generation  of  simulated  data  for  use  in  the  calculations.  The  next  step 
in  the  development  is  the  computation  of  some  examples  for  various  sampling 
schemes  and  sensor  geometries.  The  next  part  of  the  thesis  includes  a  more 
complicated  noise  model  in  the  simulation  of  the  data  and  in  the  covariance 
calculation.  Following  this  is  a  discussion  of  the  estimation  of  the  bearing  angle  using 
various  techniques  and  also  using  the  density  functions  calculated  previously.  A 
bound  on  the  estimate  of  the  variance,  the  Cramer-Rao  bound,  is  calculated  and 
compared  to  statistics  calculated  from  the  density  functions.  The  last  part  of  the 
diesis  is  a  review  of  the  results  obtained  and  recommendations  for  further  research  in 


1.2  Historical  Perspective 
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The  problem  of  determining  the  direction  of  arrival  of  wavefronts  generated  by  a 
source  in  the  presence  of  noise  has  been  discussed  extensively  in  the  literature.  The 
earliest  beginnings  lie  as  far  back  as  Leonardo  de  Vinci  who  noted  that  one  can  hear 
distant  ships  by  placing  one  end  of  an  air  filled  tube  into  the  water  and  placing  the 
other  to  the  ear. 

The  onset  of  WWI  really  began  the  explosion  of  development  of  acoustical  devices  for 
detecting  and  localizing  a  source.  The  British  and  French  are  credited  with  the  fist 
application  of  acoustical  tracking  devices  to  detect  enemy  aircraft  and  zeppelins.  The 
developments  made  in  these  areas  were  later  applied  to  the  detection  of  underwater 
sound  sources.  Early  devices  were  connected  directly  to  a  human  operators  ear  via  a 
stethoscope.  Devices  such  as  the  American  SC  tube  and  the  MB  tube  are  examples  of 
early  successful  detection  and  localization  equipment  (Burdic,  1989,  pp. 22-23). 

The  development  of  vacuum  tube  technology  during  WWII  helped  to  improve  the 
performance  of  systems  designed  for  detection  purposes.  Also  during  the  war,  the 
work  of  Noibert  Wiener  and  S.O.  Rice  in  brought  communications  theory  into  a 
new  era.  Methods  of  separating  a  signal  from  the  noise  in  which  it  is  embedded 
began  to  mature  with  the  work  of  Shannon,  Gabor,  and  Woodwind  (Burdic,  1989, 
pp.22-23). 
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With  the  end  of  WWII,  the  onset  of  the  cold  war  spawned  developments  such  as  the 
atom  bomb,  ballistic  missiles,  and  the  submarine-launched  missile.  These  new 
systems  made  the  detection  and  tracking  of  sources  a  major  priority  throughout  the 
world. 

With  the  development  of  the  digital  computer  through  the  late  60’s  and  the  70’s,  the 
use  of  the  computer  has  fully  taken  over  the  job  of  tracking.  Hence,  many  algorithms 
have  been  developed  to  solve  the  direction-of-arrival  problem.  Three  main  types  of 
algorithms  have  been  developed.  The  first  are  known  as  "signal  subspace”  methods. 
The  relationship  between  the  array  geometry  and  the  direction  of  arrival  can  be 
characterized  by  a  clearly  defined  structure  of  the  covariance  matrix.  The  signal 
subspace  methods  exploit  the  structure  of  the  covariance  to  determine  the  direction  of 
arrival.  The  second  class  of  direction  of  arrival  algorithms  are  known  as 
beamforming  techniques.  In  these  methods,  sensor  weights  and  tapped  delays  are 
used  in  order  to  alter  the  sensitivity  of  the  array  in  specific  directions.  The  way  in 
which  the  weights  are  selected  determines  the  nature  of  the  algorithm.  A  third  class 
of  algorithms  are  called  spectral  estimation  techniques.  Spectral  estimation  techniques 
utilize  the  relationship  between  estimates  of  discrete  frequencies  in  a  time  series  and 
spatial  wavenumbers  associated  with  an  interference  pattern.  These  techniques  often 
base  the  model  of  the  received  signal  on  a  state  space  model  and  adopt  the  model  to 
match  the  signal. 
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In  the  Bayes’  approach  to  estimation,  a  cost  function  which  is  a  function  of  the  error 
is  optimized  using  statistics  computed  from  the  conditional  a  posterior  probability 
function.  When  the  signal  and  noise  processes  are  Gaussian,  and  the  signal 
parameters  can  be  accurately  estimated,  the  minimum  variance  and  maximum 
likelihood  estimators  closely  approximate  each  other.  When  the  SNR  is  low  or  when 
the  array  geometry  or  medium  contain  uncertainties,  it  is  necessary  to  obtain  the 
direction  of  arrival  density  function  conditioned  on  the  actual  data  at  the  sensor 
outputs.  With  the  directly  calculated  density  function,  algorithms  can  be  used  to 
optimize  statistics  of  the  direction  of  arrival  estimate. 

1.3  Goals  of  the  Thesis 

There  are  three  main  goals  of  the  thesis.  The  first  goal  is  to  present  the  rigorous 
development  of  the  theory  involved  in  the  calculation  of  the  probability  density 
function  of  direction  of  arrival  conditioned  on  the  data  sampled  at  an  array  of  sensors. 
The  second  goal  of  the  thesis  is  to  show  how  the  density  functions  behave  when 
various  parameters  of  array  geometry  and  sampling  scheme  are  varied.  By  observing 
how  the  density  functions  change  as  a  function  of  bearing  it  is  shown  that  the  density 
functions  of  bearings  near  0°  appear  Gaussian  while  densities  near  ±90"  are  far  from 
Gaussian.  The  third  and  final  goal  of  the  thesis  is  to  give  a  detailed  discussion  of  the 
estimation  of  direction  of  arrival  using  the  densities  calculated  above. 


The  variances  of  the  estimates  of  bearing  are  then  compared  to  the  Cramer-Rao 
bound.  This  is  done  in  order  to  give  a  comparison  of  performance  of  a  Bayes’ 
optimal  estimator  with  the  C-R  bound. 


Chapter  2  Calculation  of  the  Density  Function 
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2.1  Introduction 

In  this  chapter  an  analytical  formula  is  developed  for  the  probability  density  of  the 
direction  to  a  source  in  the  far  field  of  a  array  of  sensors.  The  density  function  is 
conditioned  on  a  finite  set  of  sensor  outputs.  Both  the  signal  and  corrupting  noise  are 
assumed  to  be  zero-mean  Gaussian  random  processes.  First,  a  model  of  the  signal  is 
developed.  This  development  includes  the  sensor  geometry  and  sampling  scheme. 

The  development  of  die  signal  model  begins  with  die  signal  alone  then  is  expanded  to 
include  the  noise  process.  It  is  shown,  through  a  discussion  of  die  correlation 
function,  how  the  signal  and  noise  processes  are  characterized.  The  development  of 
the  observation  model  is  discussed  next  and  goes  through  four  stages.  First,  the 
continuous  time  signal  received  at  each  sensor  is  discussed.  This  leads  to  a  discussion 
of  the  signal  received  at  each  sensor  in  terms  of  the  signal  received  at  a  reference 
location.  The  signal  received  at  die  reference  location  is  related  to  the  signal  at  each 
of  die  season  by  the  propagation  delay  between  senson.  Third,  the  time  sampling  of 
the  signal  is  included  in  the  observation  model.  Fourth,  the  inclusion  of  additive 
noise  is  discussed.  After  die  noise  characterization  model,  the  density  function  of  the 
signal  conditioned  on  the  direction  of  arrival  is  presented.  Through  the  application  of 
Bayes'  Rule  to  this  density,  the  density  function  of  direction  of  arrival  conditioned  on 
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the  set  of  observations  is  produced.  The  most  involved  portion  of  the  development  of 
the  above  density  function  is  the  development  of  the  covariance  matrix.  A  detailed 
discussion  of  the  covariance  matrix  is  presented.  This  includes  the  relationship 
between  the  covariance  matrix  and  the  autocorrelation  of  the  signal  received  at  the 
reference  location.  Finally,  a  description  of  an  algorithm  for  the  construction  of  the 
covariance  and  conditional  density  function  is  presented. 

In  the  following  discussion  the  notation  to  be  used  is  as  follows.  The  probability 
associated  with  an  event,  A,  will  be  denoted  Prob(A).  Random  variables  will  be 
denoted  by  upper  case  letters  and  the  result  of  an  experiment  will  be  denoted  by  a 
lower  case  letter.  A  probability  density  is  denoted  by  a  lower  case  p,  with  an  upper 
case  subscript  to  denote  the  random  variable.  For  example,  px  is  the  probability 
density  associated  with  the  random  variable  X,  and  px(a)  is  the  probability  density 
function  of  the  random  variable  X  evaluated  at  a.  A  conditional  random  variable  is 
denoted  by  a  vertical  bar.  For  example,  px  )Y  is  the  probability  density  of  the  random 
variable  X  conditioned  on  the  observation  Y. 

2.2  Signal  Model 

The  purpose  of  this  section  is  to  present  the  development  of  the  covariance  matrix.  In 
order  to  fully  develop  the  covariance,  a  model  of  the  signal  and  explanation  of  the 
sensor  geometry  is  necessary.  It  is  assumed  that  there  is  a  source  producing  an 
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acoustic  signal  located  in  a  medium.  This  signal  can  be  a  complicated  combination 
of  signals  from  a  complex  source.  Examples  would  include  machinery  noise,  bearing 
or  rotation  noise,  flow  noise,  or  a  complicated  combination  of  all  of  these.  The 
source  could  produce  a  random  signal  in  Gaussian  background  noise  or  a  sinusoid 
with  high  signal  to  noise  ratio.  The  only  assumption  made  on  the  signal  and  noise  is 
that  both  are  zero  mean,  Gaussian  processes. 

The  signal  and  noise  are  received  by  an  array  of  sensors.  The  array  can  be  a  two- 
dimensional  array  of  n  sensors  or  a  linear  array.  For  die  purpose  of  simplicity,  the 
rest  of  the  discussion  assumes  a  linear  array  of  n  elements.  The  source  discussed 
above  is  assumed  to  be  in  the  far  field  of  die  array  so  that  the  signal  arrives  in  die 
form  of  plane  waves  across  the  array  as  seen  in  Figure  1. 

The  geometry  of  the  array  is  shown  in  Figure  1  as  an  n  element  array,  with 
separation  distance  d  between  sensors.  The  signal  arrives  in  the  form  of  plane  waves 
and  arrives  at  bearing  angle 


Figure  1:  Receiver  and  Signal  Geometry 

The  signal  is  a  continuous  time  process  which  is  discretized  both  spatially  and 
temporally.  The  signal  is  spatially  ‘sampled*  due  to  the  fact  that  there  is  a  finite 
distance  between  sensors.  The  signal  is  discretized  temporally  by  digitally  sampling 
the  signal.  The  sampling  interval  can  either  fixed  or  variable  length.  The 
development  of  the  covariance  which  follows  begins  with  die  signal  at  an  individual 
sensor  and  is  expanded  to  include  multiple  samples  taken  at  all  of  the  sensors.  The 
discussion  is  further  expanded  to  include  the  noise  process  in  the  following  section. 
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2.3  Development  of  the  Covariance  Matrix 

The  signal  received  at  each  of  the  sensors  is  denoted  by  the  symbol  Sa(t) 

[n- 1,2,3,...].  St(t)  denotes  the  signal  received  at  an  arbitrary  reference  sensor.  The 
signal  at  this  arbitrary  reference  sensor  is  related  to  the  signals  received  ait  the  other 
sensors  by  die  propagation  delay  between  sensors.  The  relation  between  the  signal 
received  at  the  reference  sensor  (sensor  1)  and  at  the  other  sensor  locations  assuming 
the  geometry  in  Figure  1  is  given  by  S.(t)»S,(t-rJ  where 


0) 

•  c 

d  is  the  distance  between  neighboring  sensors,  $  is  the  bearing  angle,  and  c  is  die 
sound  speed  in  the  medium.  While  only  a  linear  array  of  n  elements  is  shown  and 
planar  wavefronts  are  assumed,  the  results  presented  can  be  readily  generalized  to  a 
two  dimensional  array  with  the  sensors  having  arbitrary  separation  distance  and  the 
signal  arriving  with  curved  wavefronts  by  modification  of  (1). 

Assuming  that  die  signal  is  zero  mean  and  Gaussian,  the  process  is  characterized  by 
its  autocorrelation  function  Rs,(t,).  Given  the  autocorrelation  Rsi(t|),  pair-wise  joint 
processes  of  the  form  (S,(t,) ,  Sj(t,)]  can  be  characterized  by  the  matrix  valued 


autocorrelations: 
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Rsi(ti~xi*ti~xP 


*tfii-*f‘rx) 

XsfirW rxP 


(2) 


The  process  [Sj(t,) ,  Sj(t,)]  is  a  model  for  sampling  of  the  signal  mice  at  each  on  the 
sensors  i  and  j  at  times  t,.  Equation  (2)  is  arrived  at  by  noting  that  the  signal 
received  at  a  given  element  is  related  to  the  signal  at  the  reference  sensor  by  a  time 
shift.  The  dependance  on  can  be  incorporated  into  (2)  by  the  relation  (1).  The 
autocorrelations  (2)  can  be  extended  to  joint  processes  of  a  signal  received  at  a 
number  of  sensors  n.  Here  tt  denotes  die  time  at  which  the  sensors  are  sampled.  The 
vector  valued  process  can  be  characterized  by:  S(ti)A[S|(t|),S2(t,),...,SB(t,)]T  which 
has  the  autocorrelation  function: 


^r,(*i"Ti»*i_T»)  •••••  Rsft i~x»**i~xP 


(3) 


On  examination  of  (3)  it  is  seen  that  the  autocorrelation  is  parameterited  by  {t,}  and 
r..  Explicit  dependence  in  (3)  on  the  bearing  angle  ^  is  accomplished  via  (1). 
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The  probability  density  function  of  the  random  variable  S(tt)  conditioned  on  the 
bearing  angle  ^  is  given  by: 


It  is  seen  that  the  inverse  of  the  covariance  (3)  is  needed  in  (4).  There  are  a  number 
of  possible  bearings  for  which  the  inverse  of  the  covariance  does  not  exist.  In  the 
following  discussion,  when  a  density  function  of  the  form  (4)  appears,  its  existence 
depends  on  the  value  of  bearing  ^  and  die  value  of  the  density  at  is  understood  to 
be  defined  by  a  limit  taken  as  approaches  \ft0  since  ps(  *  fails  to  exist  only  for 
isolated  values  of  \f/.  It  is  also  recognized  that  when  a  density  function  of  the  form 
(4)  appears,  it  will  appear  in  an  expression  involving  an  integral  over  a  quadrant  of  its 
domain.  This  characterizes  a  distribution  of  the  form 

4 

Prob(X<a)  =  jpxQc)dx  (5) 

and  is  always  well  defined. 

The  covariance  (3)  must  now  be  generalized  further  to  include  multiple  samples  taken 
at  each  sensor.  The  covariance  (3)  characterizes  the  case  of  one  sample  taken  at 
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each  sensor.  For  the  purposes  of  the  work  done  here,  the  same  number  of  samples 
are  taken  at  each  sensor  and  the  samples  are  taken  at  uniform  sampling  intervals.  In 
general,  however,  the  samples  need  not  be  taken  at  the  same  time  nor  must  the  same 
number  of  samples  be  taken  at  each  sensor.  The  set  of  sampled  values  described  by 
the  vector  S  can  now  be  written: 


S(t)  -  [5,(0 

»••••  5j(0 » ••• » » ••• 


(6) 


Here  t=[t,  ,  ...tp]  is  the  set  of  p  times  that  each  sensor  is  sampled,  and  S  is  a  finite, 
jointly  distributed,  Gaussian  random  variable.  The  covariance  of  S=S(t)  is  given  by 


^j(t)-£[SSrl  (7) 

which  can  be  determined  by  using  the  autocorrelation  of  the  reference  sensor 
Rsi(t,)  at  the  appropriate  times.  The  density  function  of  S  conditioned  on  ^  is: 

**  — t — rEt*-Fr*5"<*)  ®  (» 

(2*)1  l«sl  ’ 


where  m  is  the  total  number  of  samples  taken  at  all  of  the  sensors. 
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2.4  Inclusion  of  Additive  Noise  and  the  Density  Function 

To  further  generalize  the  discussion,  the  inclusion  of  an  additive  noise  process  will 
now  be  discussed.  The  model  for  the  signal  at  the  i’th  sensor  is 

Yft)  -  Sfi)  ♦  JV/F)  (9) 

The  noise  process  N(t)*[N1(t,)...N1(g,NJ(tl)...N2(g,N,(tl)...N,(g]T  is  assumed  to  be 
a  jointly  Gaussian  process.  Since  the  signal  and  noise  are  assumed  to  be  jointly 
Gaussian,  the  joint  process  [ST(t),NT(t)]T  is  characterized  by  the  autocorrelation 
function 


*(*;♦) 


GO) 


Sampling  the  sensor  outputs  as  described  by  (9)  is  an  observation  of  the  joint  process 
described  by  (10).  Using  (6)  to  describe  the  sampled  signal  and  a  similar  expression 
to  define  the  sampled  noise  process,  the  density  of  the  finite  random  variable 
[ST,NT)T  is  characterized  by  the  large  covariance  matrix  which  includes  time  samples 


\[**m 


(id 
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In  general,  R§§ 'GA)  is  singular  for  certain  values  of  \f/.  Reasons  for  this  are  presented 
in  Chapter  4.  The  submatrix  R^M  may  contain  correlated  elements,  but  usually 
contains  an  uncorrelated  component.  ^SN  and  Rns  will  be  zero  if  the  signal  and  noise 
are  uncorrelated.  The  density  function  of  S  and  N  conditioned  on  ^  is  given  by 


1 


<2  *mir)| 


T«p  -±[StNt] 
1 


J  TR'lM 


(12) 


The  total  input  received  is  defined  as  the  random  variable 


Y-S  +  N  <13> 

The  covariance  of  Y  is  obtained  from  (11)  using  (13)  in  the  following  manner.  The 
vector  X  is  defined  as 


m 

Using  the  definition  of  X  (13)  can  be  written  as 

Y  =  [I7]X 


(14) 


(15) 


The  covariance  of  the  data  is  given  by 


16 


Kjj*£t?i'r] 

or  combining  (14)  and  (15) 


(1«) 


[/  /]  JQ 


(17) 


Where  E[XX]  is  the  covariance  of  the  vector  defined  by  (14).  Therefore,  (17)  can  be 
written  as 


(18) 


or 


^rr  M  ^ss*^ns*^sn+^nn 

Using  Bayes*  Rule  with  (19)  the  probability  density  function  for  direction  of  arrival 
conditioned  on  the  actual  samples  of  the  sensor  outputs  is  obtained.  The  result  is 


Here  p*(^)  is  a  prior  density  function  on  bearing  which  is  in  general  not  Gaussian. 
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2.5  Computation  of  the  Covariance  Matrix 

In  this  section  the  structure  of  the  covariance  R(^)  appearing  in  (12)  is  discussed.  It 
was  stated  that  die  covariance  (11)  is  constructed  by  evaluation  of  the  autocorrelation 
of  the  signal  received  at  the  reference  sensor  at  the  appropriate  times  t,  ,  tj.  In  this 
section  a  matrix  of  the  appropriate  delays  will  be  constructed  by  closely  examining 
(7).  By  selecting  values  of  the  autocorrelation  of  the  signal  received  at  the  reference 
sensor  which  correspond  to  the  delays  in  the  delay  matrix,  the  covariance  is 
constructed. 


It  was  shown  that  R(^)  *  EfSS1)  where  the  signal  S  is  given  by  (6).  Expanding  the 


matrix  notation  we  see 


s, («,»,»,)  s.d.wy 
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SST 


SJUpW 


MM 


•Ml 


•Ml 


....  ^(f,)^,) 

-  s,( rpSjCf,) 


(21) 


...  5,0^) 
..  s.(^(g 


Taking  the  expected  value  of  (21)  will  produce  the  covariance.  In  order  to  evaluate 
the  covariance  it  is  first  necessary  to  put  each  individual  element  in  (21)  in  terms  of 
the  reference  sensor  autocorrelation  Rsi(t, ,  tj).  In  order  to  better  illustrate  the  point 
of  creating  the  matrix  of  delays  a  few  examples  will  be  given.  In  order  to  simplify 
the  notation,  it  is  assumed  that  the  signal  S,(t)  is  stationary  and 
Rs,(t, ,  tj)  *  Rsi(0,t,-ta)  or  the  autocorrelation  is  if  the  form  Rm(t). 


EL  W,)]  -  EBWMhm  «  *Sl(0)  (22) 


(23) 


£[Si(W(*i)]  «  £[S,(f,  +**)$,(*,)]  -  RSi(U)  (24) 

In  the  Equations  (22),  (23),  and  (24)  the  symbol  RS1  denotes  the  autocorrelation  of  the 
signal  at  the  reference  sensor.  The  symbol  At  is  the  sampling  interval  in  seconds. 

The  symbol  r  j  i  is  the  propagation  delay  between  the  sensors  i  and  j  for  a  given 
direction  of  arrival  is  given  by 


dusm(+) 

c 


(25) 


Completing  the  above  for  all  of  the  elements  in  the  delay  matrix  and  denoting  the 
delay  matrix  D  results  in  the  matrix 
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0  Afj  2Alj  •«•••  Pi  At  i  tjj  ^i2^2Atj  .  tju+^A^ 

-Afj  0  Afj  ..... 

-2Af,  -Af,  0 . 


~pM 

~xu 

~'C12_At2 


By  matching  the  delays  in  the  above  matrix  (26)  with  the  delays  in  the  autocorrelation 
of  the  reference  sensor  Rs,(r)  the  covariance  matrix  is  constructed. 


It  is  seen  that  die  delay  matrix  is  skew  symmetric  about  a  main  diagonal  of  zeroes. 
The  covariance,  however,  will  be  symmetric  due  to  the  symmetry  of  the 
autocorrelation.  The  covariance  R(^)  is  obtained  by  addition  of  to  R&  assuming 
that  the  signal  and  noise  are  uncorrelated  so  that  R&  and  R*§  are  zero.  More 
discussion  of  the  covariance  is  given  later  in  relation  to  specific  noise  models. 
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2.6  Description  of  the  Code 

In  order  to  calculate  the  density  (20),  a  computer  program  was  written  to  compute, 
for  each  the  covariance  RyT(^),  the  inverse  RyT  ,(^),  and  the  determinant 
DetCRyT(^))  and  pre-  and  post-multiply  die  inverse  by  the  vector  of  sampled  values. 
The  possible  bearings  are  considered  for  a  linear  array  of  sensors  with  equal  spacing 
which  are  sampled  at  die  same  time  with  die  same  sampling  intervals  at  each  sensor. 
The  range  of  possible  bearings  is  -t/2  to  x/2  radians  since  for  a  linear  array, 
reciprocal  bearings  are  expected  and  a  "forward  looking”  array  will  be  considered. 

The  covariance  and  its  inverse  must  be  constructed  in  order  to  calculate  die  density 
(20).  In  order  to  accomplish  this,  two  main  steps  are  taken.  The  first  is  to  construct 
a  matrix  of  delays  which  are  based  on  the  propagation  delays  between  sensors  and  the 
sampling  interval  at  die  sensor  outputs.  The  second  is  to  match  the  delays  in  the 
delay  matrix  to  die  delays  in  the  autocorrelation  of  the  reference  sensor  and  place  the 
corresponding  autocorrelation  in  the  position  of  die  matching  delay. 

The  matrix  of  delays  is  constructed  from  groups  of  submatrices.  Symmetry  within 
die  covariance  greatly  reduces  the  number  of  actual  elements  in  the  delay  matrix  to  be 
filled.  For  the  case  of  equal  numbers  of  samples  taken  at  all  of  die  sensors,  the 
covariance  matrix  is  a  square  matrix  with  dimensions  [n*P  x  n*P]  where  n  is  the 
number  of  sensors  and  Pis  the  number  of  samples  taken  at  each  sensor.  Upon  close 
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examination  of  (26)  it  is  seen  that  the  subblocks  which  make  up  the  main  diagonal  are 
identical.  The  subblocks  containing  the  main  diagonal  are  calculated  from  the 
sampling  interval  and  number  of  samples  taken  at  the  reference  sensor.  This  means 
that  there  are  no  cross-terms  in  the  main  diagonal  subblocks  and  that  the  delay 
between  the  sensors  given  by  (25)  does  not  appear  in  the  main-diagonal  subblocks.  It 
is  necessary  only  to  calculate  one  such  subblock  and  to  place  it  into  positions  along 
the  main  diagonal  of  the  delay  matrix.  In  calculating  the  main-diagonal  subblocks,  it 
is  necessary  due  to  symmetry  only  to  fill  the  elements  which  lie  below  the  main 
diagonal  of  zeroes.  By  filling  the  elements  below  the  main  diagonal  with  the  proper 
delay  values  and  all  other  elements  with  zeroes,  the  subblock  is  calculated  by  adding 
the  negative  of  the  transpose  of  die  half-filled  subblock  to  the  half-filled  subblock.  It 
should  be  noted  the  number  of  subblocks  is  n*  where  n  is  again  the  number  of 
sensors.  It  is  also  seen  that  die  subblocks  which  lie  above  the  main  diagonal  of 
subblocks  are  the  transposes  of  those  lying  below  the  main  diagonal  subblocks. 
Equation  (27)  gives  an  example  of  the  delay  matrix  for  two  sensors  in  a  linear  array 
taking  four  samples  at  a  rate  of  t,. 
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TU  *12-',  ti2+2*4  Ti2+3d 
Xi2~(»  XI2  *12**s  TU+^ 
Xi2-^fl  X12"f»  TU  T12*fi 
tl2"^»  ti2~2f,  tu”f»  tu  (27) 
0  r,  2la  3fJ 
**,  0  f,  2X, 

|  '*#  'fi  0  li 

[“3 1,  -2 1,  -t,  Oj 

Once  die  delay  matrix  is  filled,  die  covariance  is  constructed  by  matching  the  delays 
in  the  delay  matrix  to  the  delays  in  the  autocorrelation  of  die  signal  received  at  the 
reference  sensor.  The  autocorrelation  of  the  signal  at  the  reference  sensor  is  sampled 
and  stored  in  a  data  vector.  In  practice,  the  matching  is  accomplished  by  executing  a 
binary  search  mi  the  data  vector  for  each  delay  put  into  the  delay  matrix  (27).  This 
is  done  so  that  die  symmetry  of  the  autocorrelation  can  be  exploited  to  reduce  the 
computation  time,  since  only  half  of  the  elements  must  be  searched. 

The  above  must  be  completed  for  each  bearing  angle  of  interest  since  the  propagation 
delay  between  sensors  used  in  constructing  the  covariance  is  a  function  of  bearing 
angle. 


2.7  Summary  of  Chapter  2 
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In  this  chapter  the  formulation  of  the  probability  density  function  of  direction-of- 
arrival  conditioned  on  the  actual  sample  values  taken  at  each  of  the  sensors  has  been 
developed.  The  result  was  arrived  at  by  development  of  the  density  function  of  the 
signal  conditioned  on  die  bearing  angle  and  the  use  of  Bayes*  rule  to  arrive  at  the 
density  function  of  bearing  angle  conditioned  on  the  sampled  signal.  The 
development  of  the  density  function  began  with  the  density  function  conditioned  on 
the  signal  and  was  expanded  to  include  the  case  of  numerous  samples  taken  at  each 
sensor.  This  was  finally  expanded  to  include  the  case  of  additive  noise  in  the  output 
of  the  sensors.  The  case  of  a  linear  array  of  sensors  was  presented,  with  the 
possibility  of  expanding  die  result  to  include  two  dimensional  arrays  of  sensors. 

Following  the  development  of  the  density  function  a  discussion  of  how  the  covariance 
is  actually  constructed  is  presented.  This  lead  to  a  presentation  of  the  algorithm  used 
to  calculate  the  density  function  given  an  observation  vector  which  contains  sample 
values  of  the  signal  at  all  of  the  sensors  in  the  array. 


Chapter  3  Discussion  of  Simulation  and  Examples  of  Density  Functions 


3.1  Introduction 

In  this  chapter,  several  illustrative  examples  of  probability  density  functions  of 
bearing  angle  conditioned  on  the  sampled  signal  data  are  explored.  The  general 
model  for  an  array  of  sensors  is  given  in  Figure  2. 


Figure  2:  Model  for  Sensor  Array 
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In  Figure  2  a  linear  array  of  n  sensors  is  shown.  The  block  labeled  H(u)  represents  a 
linear  filter  with  frequency  response  H(«).  This  filter  can  be  a  high-pass,  low-pass, 
or  band-pass  filter.  The  N/s  represent  input  noises  which  are  uncorrelated  between 
sensors.  The  N0’s  represent  output  noises  which  are  also  uncorrelated  between 
sensors.  The  examples  presented  in  the  following  sections  use  variations  on  the 
sensor  model  presented  in  Figure  2. 

For  the  examples  given,  the  signal  used  is  a  narrow-band  signal  produced  by  a 
MATRIXxtm  simulation  shown  in  Figure  3.  The  prior  density  function  as  given  in 
Equation  (18)  on  bearing  is  assumed  to  be  uniform  over  the  range  -90°  to  90”.  The 
code  which  generated  the  following  plots  is  given  in  Appendix  A.  First,  a  discussion 
of  the  simulation  of  the  signal  and  noise  at  each  of  the  sensors  is  given.  This  will 
include  a  presentation  of  the  properties  of  a  narrow  band  filter  which  is  used  in  the 
sensing  scheme.  Following  a  discussion  of  the  signal  and  noise  simulation,  a  short 
discussion  of  the  definition  of  signal-to-noise  ratio  is  given. 

Examples  to  be  discussed  will  include  the  following.  The  first  example  is  the  case  of 
an  array  of  two  sensors  with  one  source  located  at  different  bearings  in  the  half  space 
in  front  of  the  array.  The  second  example  is  an  array  of  two  sensors  and  one  source 
where  the  spacing  of  the  sensors  is  varied  as  a  function  of  the  source’s  center- 
frequency  wavelength.  The  third  example  is  for  an  array  of  two  sensors  and  one 
source  where  the  number  of  samples  taken  per  sensor  is  varied.  The  fourth  example 
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is  for  an  array  of  two  sensors  where  the  signal  to  noise  ratio  is  varied.  For  all  of  the 
examples  except  for  the  second,  the  sensor  separation  distance  is  one-fourth  of  the 
source’s  center-frequency  wavelength. 


3.2  Simulation  of  the  Signal  and  Noise 

It  has  been  shown  in  the  previous  sections  that  the  density  function  for  the  direction 
of  arrival  conditioned  on  the  actual  sensor  outputs  can  be  determined.  In  this  section, 
the  method  employed  to  generate  the  data  which  is  to  be  sampled  is  discussed. 

A  commercially  available  software  package,  MATRIXx™  is  employed  to  do  the 
simulation.  A  linear  system  is  driven  with  a  Gaussian,  uncorrelated  sequence  to 
simulate  broadband  noise.  Further  discussion  of  the  simulation  outputs  is  given  in 
Section  3.4.  This  produces  a  narrow-band  signal  with  slowly  varying  phase.  For  the 
examples  in  this  section,  only  one  linear  system  is  used  for  simplicity.  Many  such 
systems  can  be  used,  however,  and  their  outputs  summed  to  produce  a  more 
complicated  signal.  For  example  a  signal  with  wide-band  characteristics  can  be  used 
to  simulate  thermal  or  flow  noise  at  the  sensor  inputs.  A  narrow-band  signal  or  many 
narrow-band  signals  may  also  be  present  in  the  signal.  These  could  be  models  for 
rotating  machinery  that  is  slightly  out  of  balance.  For  the  more  simple  signal  used  in 
the  following  examples,  the  frequency  response  of  the  filter  can  be  seen  in  Figure  3. 
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The  output  of  the  narrow-band  filter  is  split  into  different  channels  which  represent 
the  different  sensors.  For  all  of  the  examples  discussed  in  this  chapter,  the  reference 
sensor  is  represented  by  channel  1.  The  signals  which  go  into  the  other  channels  are 
delayed  by  an  amount  given  by  Equation  (1)  to  simulate  the  propagation  delay 
between  sensors.  To  the  output  of  each  channel,  Gaussian  noise  is  added.  The  noise 
at  each  sensor  is  made  to  be  uncorrelated  with  the  noise  at  other  sensors  by  changing 
the  seed  in  the  random  noise  generators.  A  schematic  for  the  signal  and  noise 
generators  is  given  in  Figure  3  for  the  case  of  one  source  and  two  sensors. 


whit*  Nois* 
Generator 


Figure  3:  Schematic  for  Simulation  of  Data 
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The  Gains  G  shown  in  Figure  3  are  installed  as  a  means  of  controlling  the  SNR.  The 
gains  are  further  discussed  in  the  next  section.  The  simulation  is  run  for  a  time  long 
enough  for  the  signal  to  be  uncorrelated  with  the  initial  state.  The  simulation  in 
Figure  3  produces  three  vectors  of  data.  The  three  data  vectors  are  the  signal 
received  at  the  two  sensors  given  by  outputs  1  and  2  and  the  autocorrelation  of  the 
signal  at  the  reference  sensor  which  is  computed  after  running  the  simulation  by 
correlating  the  output  3  with  itself.  This  is  the  autocorrelation  of  the  signal  at  the 
reference  sensor  discussed  in  Section  2.3. 

The  outputs  of  the  two  channels  are  then  sampled  at  a  given  sampling  interval.  Since 
the  simulation  produces  discrete  data,  the  data  is  sampled  at  a  rate  which  is  a  multiple 
of  the  output  interval.  The  sampled  values  are  ordered  into  two  data  vectors  which 
represent  the  data  received  at  the  two  sensors.  These  data  vectors  are  then  combined 
into  the  signal  vector  given  an  expression  similar  to  that  of  (6)  except  that  the  noise  is 
now  added.  A  sample  of  the  output  of  the  simulation  from  channel  1  is  seen  in 
Figure  4.  The  data  created  by  the  simulation  are  the  signal,  as  sampled  at  the 
sensors,  and  the  autocorrelation  of  the  signal  at  the  reference  sensor.  This  data  is 
then  passed  to  a  FORTRAN  subprogram  to  calculate  the  covariance  matrix  discussed 
in  Section  2.3  and  the  probability  density  conditioned  on  the  sensor  outputs. 
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Figure  4:  Output  of  MATRIXx  Simulation  for  SNR =5. 
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Figure  5:  Frequency  Response  of  the  Filter  referenced  to  rad/sec. 


3.3  Signal-to-Noise  Ratio 
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In  this  section  the  definition  of  signal-to-noise  ratio  is  given  in  relation  to  the  signal 
model  discussed  in  the  previous  section.  In  particular  the  relation  between  signal-to- 
noise  ratio  and  the  autocorrelation  of  the  signal  and  noise  components  of  the  output  at 
channels  1  and  2. 


The  signal-to-noise  ratio  is  defined  as  the  ratio  of  the  signal  and  noise  variances. 


sm.  Signal  Variance 
Noise  Variance 

The  signal  variance  (as2  )  in  (27)  becomes 


07) 


=  /  b-pjfpis)  ds  *  t,2-p,2  08) 

or 
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and 

o.* 1*.’  (3#) 


For  the  purpose  of  this  work  the  signal  and  noise  are  assumed  to  be  zero-mean  and 


Gaussian.  The  signal-to-noise  ratio  is  therefore  given  by 
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SNR 
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(31) 


3.4  Analysis  of  Simulation  Outputs 

In  this  section  an  analysis  of  the  outputs  on  the  MATRIXx™  simulation  is  presented. 
The  ’whiteness’  of  the  noise  generator  driving  the  linear  system  is  at  best 
questionable  and  criteria  must  be  established  to  ensure  that  the  noise  generator 
appears  white  to  the  linear  system.  Another  area  to  be  investigated  is  how  the 
signal-to-noise  ratio  as  presented  in  the  previous  section  is  set  in  the  simulation.  The 
discussion  begins  with  an  analysis  of  the  white  noise  generator. 

The  noise  generator  which  drives  the  linear  system  in  Figure  3  is  in  reality  a  Gaussian 
random  number  generator.  The  user  sets  the  standard  deviation,  the  seed,  and  the 
output  interval.  The  output  of  the  generator  typically  looks  like  the  signal  shown  in 
Figure  6. 


The  autocorrelation  (32)  looks  like 
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Figure  7:  Autocorrelation  of  the  Noise  Sequence. 

Taking  the  Fourier  Transform  of  the  autocorrelation  R^r)  the  autospectrum 
is  obtained  (Bendat  and  Piersol,  1986,  p.  121) 

m 

V.(u)  ‘  /  *T°»WoT,'")  <33) 

Hi  “  ^  * 

This  function  has  a  main  lobe  the  width  of  which  is  controlled  by  the  output  interval 
of  the  noise  generator  At  It  is  seen  that  by  suitably  choosing  the  value  of  At  the 
main  lobe  can  be  made  wide  enough  so  that  the  spectrum  at  the  input  to  die  filter 


appears  white  over  some  range  containing  the  center-frequency  of  the  filter.  The 
criteria  set  for  the  noise  signal  at  the  input  to  the  filter  to  be  white  is  that  the  slope 
and  magnitude  of  the  autospectrum  remain  within  10%  of  the  center-frequency  level 
over  twice  the  filter  bandwidth.  The  filter  band-width  is  defined  as  the  frequency 
range  covering  the  -3dB  points  on  the  frequency  response  curve. 

The  frequency  response  of  the  filter  is  given  by  the  function  H(s) 


H(a) - iiii -  (34) 

1 

where  0  is  set  by  the  user.  To  see  the  response  of  the  filter,  the  function  |  H(u)  j 2  is 
plotted  vs.  <■>  in  Figure  5  for  a  value  of  0  set  at  1. 

Given  the  functions  ]  H(cj)  | 2  and  the  power  spectrum  at  the  output  of  the 

filter  is  produced  by  the  relation 


SJa)  -  |ff(«)|2S^(«) 


(35) 


and  by  taking  the  inverse  Fourier  Transform  of  Sa(«)  the  autocorrelation  of  the 
output  of  the  filter  is  obtained. 
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RJx)  «  f  SJa)  e  -J"xdu>  (36) 

The  autocorrelation  has  a  maximum  at  a  delay  value  r*0.  By  finding  the  maximum, 
the  gain  G  in  Figure  2  can  be  set  by  choosing  G  as  the  inverse  of  the  square  root  of 
die  maximum  of  R^r).  This  will  normalize  the  function  R^r)  such  that  the 
maximum  is  set  to  1. 

The  signal-to  noise  ratio  is  set  by  adjusting  the  standard  deviation  on  the  additive 
noise  generators.  Since  die  SNR  is  defined  as  the  ratio  of  die  autocorrelations  of  die 
signal  and  noise,  respectively,  and  since  die  autocorrelation  of  the  signal  at  zero  delay 
is  one,  and  die  value  of  die  noise  autocorrelation  at  zero  delay  is  the  standard 
deviation  of  die  noise,  the  SNR  is  given  by 


SNR 


JWQ) 

JWO) 


(37) 


Therefore  by  a  simple  algebraic  manipulation  of  (37)  die  SNR  is  chosen  and  the 
standard  deviation  of  die  noise  generators  are  set  accordingly. 
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3.5  Example  1:  Effect  of  Change  in  Bearing  on  the  Density  Function 

In  die  following  example,  the  geometry  of  die  array  and  source  are  as  given  in  Figure 
(1).  The  data  is  produced  using  the  simulation  discussed  in  Section  3.2  and  the 
probability  density  functions  are  calculated  by  the  code  given  in  Appendix  A.  The 
parameters  which  remain  constant  throughout  the  example  are  die  number  of  samples 
per  sensor  P.-20;  the  sampling  interval  t,*l  sec;  die  sensor  separation  distance 
d«X/4  where  X  is  die  wavelength  associated  with  the  signal  center  frequency.  A 
parameter  which  will  be  changed  is  the  actual  bearing  angle.  This  is  accomplished  by 
manipulation  of  the  delay  seen  in  the  second  channel  in  Figure  3  which  is  related  to 
die  actual  bearing  angle  by  Equation  (1).  A  plot  of  the  probability  density  function 
conditioned  on  a  set  of  sensor  outputs  for  the  above  parameters,  for  bearings  of  (P, 
20*, 45°,  75*  ,  and  85*  is  given  in  Figure  (8).  The  covariances  for  each  bearing  angle 
for  which  die  density  is  calculated  are  calculated  off-line.  The  covariance  is 
constructed  using  the  autocorrelation  of  the  signal  at  the  reference  sensor.  The  noise 
covariance  is  added  to  the  signal  covariance  in  die  form  of  a  diagonal  matrix  using 
(18).  The  main  diagonal  of  the  noise  covariance  contains  the  only  non-zero  values  in 
the  matrix.  This  occurs  due  to  the  fact  that  the  noise  is  uncorrelated  between  sensors 


and  with  itself  in  time. 
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Figure  8:  Probability  Density  Functions  for  Varying  Bearing  Angles 
SNR  =5,  10  Samples  per  Sensor  at  2  sec  Intervals. 


The  density  functions  seen  in  Figure  8  are  typical  examples  for  each  given  bearing. 
The  general  trends,  when  large  numbers  of  density  functions  are  generated  for 
different  data  sets  are  shown.  It  is  possible,  however  that  for  a  given  bearing,  the 
density  function  can  look  quite  different  due  to  locally  higher  noise  levels  or  a  variety 
of  other  parameters.  An  example  of  varying  density  functions  for  different  data  sets 
are  given  in  Figure  9.  The  densities  are  slightly  different  for  each  data  set.  This 
may  occur  due  to  the  fact  that  small  data  sets  are  used  and  may  contain  locally  higher 


or  lower  noise  levels. 
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Bearing  Angle  <**> 


Figure  9:  Varying  Density  Functions  For  Different  Data  Sets.  SNR =5 
10  Samples  per  Sensor  at  2  Sec  Intervals 

The  density  function  for  zero  bearing  appears  Gaussian.  Typically,  as  the  bearing  of 
the  source  approaches  ±90°,  the  density  functions  tend  to  be  broader  with  lower 
peaks.  This  indicates  that  less  accurate  estimates  can  be  made  on  the  bearing  angle  as 
the  bearing  approaches  ±90°.  Very  near  to  ±90°  the  density  functions  begin  to 
narrow  and  the  peaks  tend  to  get  higher.  This  occurs  since  the  density  functions  are 
essentially  "windowed"  by  the  prior  density  in  (18).  The  curves  also  tend  to  be  more 
skewed  as  the  bearing  of  the  source  approaches  ±90°  indicating  that  the  disparity 
between  the  means  and  maximums  of  the  densities  is  growing  larger. 
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3.6  Example  2:  Effect  of  Sensor  Separation  on  the  Density  Function 

In  this  section  the  effect  of  changing  the  sensor  separation  on  the  conditional  density 
function  is  discussed.  In  Section  1.2  it  was  shown  that  the  propagation  delay  between 
sensors  was  a  function  of  separation  distance.  This  relation  was  given  by  equation 
(1).  By  noting  that  the  sound  speed  in  the  medium  is  given  by  c=Xv=2tX/u  which 
gives  the  new  form  to  (1) 

r.2TrfsinQfr)  (3g) 

Xo> 

In  this  example  the  distances  d  are  chosen  as  a  function  of  X,  where  X  is  the 
wavelength  associated  with  the  center-frequency  of  the  narrow-band  signal  discussed 
in  section  3.1  and  v  is  the  center-frequency  in  Hz..  The  distances  chosen  to  illustrate 
this  example  are  d=X/4,  X/2,  X,  and  2X.  Other  parameters  are  the  number  of 
samples  per  sensor,  20,  fo=l  rad/sec,  the  sampling  interval  is  1  sec,  and  the  SNR 
is  5.  Plots  of  probability  density  functions  vs.  bearing  angle  for  the  above  separation 
distances  and  parameters  can  be  seen  in  Figures  10  and  1 1.  It  can  be  seen  from 
these  plots  that  as  separation  distance  goes  from  X/4  to  X/2,  the  density  function 
becomes  less  broad  and  has  a  higher  peak.  This  occurs  since  the  aperture  is 
effectively  widened  in  the  spatial  domain.  As  the  separation  distance  is  widened 
beyond  X/2  (referred  to  as  the  Nyquist  Separation  from  this  point  on)  it  is  seen  that 
multiple  peaks  occur  in  the  density.  This  then  appears  to  predict  a  certain  amount  of 
spatial  aliasing. 
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Figure  10:  Probability  Density  vs.  Bearing  Angle  For  Varying  Separation  Distances 
SNR-5,  10  Samples  per  Sensor  at  2  sec  Intervals. 
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Figure  11:  Probability  Density  vs.  Bearing  Angle  For  Varying  Separation  Distances 
SNR-5,  10  Samples  per  Sensor  at  2  sec  Intervals. 
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3.7  Effect  Of  The  Number  of  Samples  per  Sensor  on  the  Density  Function 

In  this  section,  the  effect  of  changing  the  number  of  samples  per  sensor  is  illustrated. 
The  number  of  samples  per  sensor  will  be  S,  10,  IS,  and  20.  The  actual  bearing  of 
the  source  is  0°,  the  separation  distance  is  X/4,  the  sampling  interval  is  1  sec,  and  the 
SNR  is  S.  Typical  density  functions  for  the  above  parameters  are  seen  in  Figure  12. 
Again,  the  density  functions  shown  reflect  the  general  trend  for  a  large  number  of 
density  functions  generated  for  different  data  sets  for  a  given  number  of  samples  per 
sensor.  It  is  possible  for  the  density  functions  to  have  different  shapes  for  the  same 
number  of  samples  per  sensor.  As  was  the  case  in 
presented  in  Figure  9. 


Figure  12:  Density  Functions  For  Varying  Numbers  of  Samples  per  Sensor 
SNR-5.  Sampling  interval- 1  sec. 
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As  can  be  seen  in  Figure  12,  the  general  trend  is  that  as  more  samples  per  sensor  are 
taken,  the  density  function  becomes  more  narrow  and  a  better  estimate  of  the  actual 
bearing  can  be  made.  The  same  is  true  for  bearings  other  than  zero  bearing,  however 
zero  bearing  best  illustrates  the  effect  of  changing  the  number  of  samples  per  sensor 
on  the  density  function. 

3.8  Effect  of  the  Signal-To-Noise  Ratio  on  the  Density  Function 

In  this  section  the  effect  of  changing  the  signal-to-noise  ratio  on  the  density  function 
is  discussed.  The  signal  used  is  the  same  as  that  used  in  the  previous  sections.  In 
this  example,  the  signal-to-noise  ratio  is  varied  through  the  values  of  SNR=2,5,10, 
and  20.  The  actual  bearing  of  the  source  is  0°.  The  separation  distance  is  X/4  where 
X  is  the  wavelength  associated  with  the  center  frequency  of  the  narrow-band  signal. 
The  number  of  samples  taken  is  20  at  a  rate  of  1  sample  per  second.  The  density 
functions  for  the  above  parameters  with  signal-to-noise  ratios  of  2,5,10,  and  20  are 
seen  in  Figure  13. 

The  density  functions  shown  illustrate  the  general  trend  for  changes  in  the  signal-to- 
noise  ratio.  It  is  possible  that  the  density  functions  will  take  on  a  different  shape  for 
other  data  sets  with  all  other  parameters  kept  constant  similar  to  the  case  as  given  in 
Figure  9. 
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Figure  13:  Density  Functions  for  a  Change  in  Signal-to-Noise  Ratio 
10  Samples  per  Sensor  at  2  sec  Intervals. 


As  is  shown  by  Figure  13,  the  density  functions  become  increasingly  narrow  as 


signal-to-ftoise  ratio  increases. 

[ 
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3.9  Summary  of  Results 

Four  illustrative  examples  have  been  explored  in  this  chapter.  The  effects  of 
changing  the  source  bearing  angle,  the  sensor  separation  distance,  the  number  of 
samples  taken  at  each  sensor,  and  the  signal-to-noise  ratio  have  been  presented.  The 
results  of  these  examples  are  summarized  here. 

It  has  been  shown  that  as  the  bearing  angle  approaches  ±90°  the  peaks  of  the  density 
functions  become  lower  and  the  curves  become  broader,  making  an  estimate  of 
direction  of  arrival  less  certain.  Beyond  this  point,  the  peaks  begin  to  narrow  again 
and  the  curves  become  more  narrow,  however  the  peaks  of  the  curves  are  less  well 
defined.  The  peaks  of  the  density  functions  become  more  flat,  thus  making  an 
accurate  estimate  of  bearing  angle  less  certain. 

By  changing  the  sensor  separation  distance,  it  is  shown  that  for  a  separation  distance 
of  X/4,  where  X  is  the  wavelength  associated  with  the  center-frequency  of  the  narrow- 
band  signal,  a  single-peak  well  defined  density  function  exists.  When  the  separation 
distance  is  expanded  to  X/2,  the  peak  becomes  higher  and  the  density  function  curve 
becomes  more  narrow.  This  signifies  that  a  better  estimate  of  the  bearing  angle  can 
be  made  with  a  wider  separation  distance  of  X/2.  This  however  is  the  widest  possible 
sensor  separation  distance  since  beyond  this  distance  multiple  possible  bearings  exist. 
This  signifies  a  kind  of  spatial  aliasing  analogous  to  aliasing  experiences  in  digital 
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sampling.  The  separation  distance  of  X/2  is  analogous  to  the  Nyquist  rate  in  digital 
signal  processing. 

For  the  case  of  changing  the  number  of  samples  taken  at  each  sensor,  it  is  shown  that 
as  the  number  of  samples  increases,  the  density  function  becomes  more  narrow  with  a 
higher  peak.  As  a  result,  a  better  estimate  of  the  direction  of  arrival  can  be  made. 

By  changing  the  signal-to-noise  ratio  it  is  shown  that  as  the  signal  becomes  more 
discemable  from  the  noise,  the  density  function  becomes  more  narrow  with  a  higher 
peak. 
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Chapter  4:  The  Conditional  Density  Function  and  Bearing  Estimates 


4.1  Introduction 

In  this  chapter,  the  main  focus  is  on  the  estimation  of  the  bearing  angle  and  other 
statistical  properties  derived  from  the  conditional  density  functions.  The  probability 
density  functions  which  are  to  be  used  in  this  chapter  are  calculated  using  the  same 
methods  as  those  in  Chapter  3.  The  simulation  which  models  the  signal  and  array  of 
sensors  will  be  different  from  that  in  Chapter  3.  In  the  next  section,  the  simulation 
which  is  used  will  be  discussed  in  detail.  In  Section  4.3,  an  analysis  of  the 
covariance  is  again  presented.  This  is  necessary  since  a  new  noise  scheme  is  used  in 
the  simulation.  Also  in  Section  4.3  is  a  discussion  of  some  of  the  limitations  involved 
in  the  methods  used  to  calculate  the  conditional  density  functions.  Following  the 
discussion  of  the  covariance,  examples  of  density  functions  calculated  for  the  new 
simulation  are  presented.  In  the  next  section,  the  examples  mentioned  above  are  used 
to  complete  a  parametric  study.  In  this  study,  such  statistical  properties  as  the 
maximum  a  posteriori  estimate,  the  minimum  mean  square  error  estimate,  and  the 
variances  on  the  estimates  are  presented.  Also  for  a  large  number  of  densities,  the 
means  and  variances  of  these  estimates  are  presented  with  confidence  intervals  on 
the  calculations.  These  estimate  variances  are  compared  with  the  Cramer-Rao  lower 


bound  for  estimate  variances. 
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4.2  Discussion  of  the  Simulation 

In  this  section,  a  discussion  of  the  MATRJXx  simulation  which  produces  the  data 
used  in  this  chapter  is  presented.  The  simulation  produces  data  for  a  two  sensor 
array.  The  data  produced  by  the  simulation  are  the  signal  received  at  each  sensor 
and  the  noise  at  each  sensor.  The  signal  is  interpreted  as  a  broadband  random 
process  which  passes  through  the  filter  and  becomes  narrow-band.  It  is  assumed  that 
the  noise  received  at  each  of  the  sensors  is  uncorrelated  with  that  received  at  the  other 
sensor.  A  schematic  of  the  simulation  is  given  in  Figure  14. 


Figure  14:  Schematic  of  Data  Simulation 
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All  of  the  noise  sources  are  uncorrelated  with  each  other.  This  is  accomplished  by 
setting  all  of  the  seeds  in  the  noise  generators  to  be  different  for  each  run  of  the 
simulation.  Hie  narrow-band  filters  seen  in  Figure  14  are  the  same  as  those  in 
Chapter  3.  The  main  difference  between  the  simulation  in  Chapter  3  (Figure  3)  and 
the  simulation  for  this  chapter  is  that  the  simulation  as  given  in  Figure  14  has  narrow- 
band  noise  that  is  uncorrelated  at  each  of  the  sensors.  A  physical  example  of  this 
would  be  flow  noise  which  is  passed  through  a  narrow  band  filter.  As  in  the 
simulation  in  Chapter  3,  the  noise  spectrum  is  sufficiently  broad  compared  to  the 
band-width  of  the  filter  that  it  can  be  assumed  white.  The  noise  sources  which  appear 
after  the  filter  are  also  broad-band  random  noise  sources  which  can  be  interpreted 
physically  as  electrical  noise.  In  all  of  the  simulation  runs,  the  electrical  noise  is  set 
to  be  considerably  lower  than  that  of  the  noise  found  external  to  the  sensor.  This 
assures  that  the  covariance  will  not  be  singular  by  adding  a  diagonal  component 

In  the  construction  of  the  covariance  the  narrow-band  noises  appear  as  block  diagonal 
elements.  There  is  no  correlation  of  the  noises  between  sensors.  This  implies  that 
the  inter-element  cross-terms  in  the  covariance  will  all  be  zero.  The  autocorrelation 
of  the  noise  at  each  sensor  has  a  non-zero  value.  Since  the  broad-band  noise  is 
modeled  as  white,  the  electrical  noise  contributes  only  to  the  main  diagonal  of  the 
covariance  matrix  (i.e.  each  sample  is  only  correlated  with  itself).  Further 
discussion  of  the  covariance  matrix  appears  in  the  next  section. 
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There  now  exist  two  noise  sources.  The  signal-to-noise  ratio  will  now  be  defined  as 
the  ratio  of  the  sum  of  the  variances  of  the  noises  to  the  variance  of  the  signal. 

4.3  Discussion  of  the  Covariance  Matrix 

In  this  section,  the  changes  which  occur  in  the  covariance  matrix  are  discussed.  As 
mentioned  above  the  main  difference  is  that  there  are  now  block  diagonal  elements  in 
the  covariance  due  to  the  filtered  noise.  Since  the  signal  and  noise  both  pass  through 
die  same  filter,  and  since  both  are  broad-band  processes,  the  autocorrelations  of  the 
signal  and  noise  will  have  autocorrelation  functions  which  have  the  same  shape,  but 
have  different  levels.  That  is,  the  noise  autocorrelation  is  a  multiple  of  the  signal 
autocorrelation.  Since  the  noise  processes  are  uncorrelated  between  sensors,  the 
noises  show  up  in  the  covariance  as  a  block  diagonal  matrix  added  to  the  covariance 
of  the  signal  alone.  That  is,  the  noise  covariance  which  appears  in  (18)  is  a  block 
diagonal  matrix  where  the  subblocks  which  appear  above  and  below  the  main  diagonal 
contain  only  zeroes.  The  main  diagonal  of  the  covariance  matrix  is  the  sum  of  the 
value  of  the  autocorrelation  of  the  signal  at  zero,  the  variance  of  the  broad-band 
noise,  and  the  variance  of  the  narrow-band  noise. 

Another  important  point  to  make  about  the  covariance  matrix  is  that  near  singularities 
occur  for  various  bearing  angles.  The  number  and  frequency  of  near  singularities  can 
be  related  directly  to  the  number  of  samples  taken  at  each  sensor,  the  sampling 


interval,  and  the  spacing  between  sensors.  The  near  sing  ’> cities  occur  for  bearings 
which  have  a  propagation  time  delay  between  sensors  (r)  which  is  equal  to  the 
sampling  interval  and  integer  multiples  of  the  sampling  interval.  Figure  IS  illustrates 
this  urng  one  period  of  a  cosine  wave  as  the  autocorrelation.  The  symbols  ,  x, 
represent  the  values  of  the  autocorrelation  of  the  signal  at  the  reference  sensor  with 
later  samples  at  the  reference  sensor.  The  symbols,  □,  represent  the  correlation 
between  the  samples  at  the  reference  sensor  and  the  samples  at  a  second  sensor.  As 
dt,  the  sampling  interval,  approaches  the  magnitude  of  the  propagation  delay,  r,  the 
near  singularities  occur  (i.e.  when  the  x’s  and  the  CPs  overlap).  This  happens  since 
the  elements  in  die  covariance  off-diagonal  subblocks  are  shifted  versions  of  the 
diagonal  subblocks.  When  this  occurs,  the  covariance  in  near  singular.  It  has  been 
shown  that  the  least  number  of  near-singularities  occur  when  the  sampling  interval 
corresponds  to  the  Nyquist  sampling  interval  corresponding  to  the  center  frequency  of 
the  signal.  Also  the  spacing  should  be  near  the  spacing  corresponding  to  a  half 
wavelength  of  the  center  frequency  of  die  signal. 
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Figure  15:  Sampling  Interval,  Propagation  Delay,  and  Covariance  Singularities 

In  order  to  illustrate  the  effect  of  covariance  singularity,  the  condition  number  of  the 
matrix  is  plotted  as  a  function  of  bearing  angle.  The  condition  number  is  defined  as 
the  ratio  of  die  largest  to  the  smallest  eigenvalues  of  the  covariance.  The  larger  the 
condition  number,  the  closer  to  singular  the  covariance  matrix.  Figure  16  shows  the 
condition  number  as  a  function  of  bearing. 

Figure  16  shows  plots  of  the  condition  number  as  a  function  of  bearing  angle  for 
three  sampling  schemes.  The  lowest  curve  is  for  Nyquist  sampling  which  is  10 
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samples/sensor  at  2  sec  intervals.  The  middle  curve  is  the  condition  number  plot  for 
20  samples  per  sensor  at  a  1  sec  sampling  interval.  The  bearing  which  corresponds  to 
the  1  sec  sampling  interval  is  ±39.5°(this  is  found  using  Equation  1).  It  is  seen  that 
the  peak  does  in  fact  occur  at  -39.5°.  The  highest  curve  is  for  over-sampling  with  20 
samples/sensor  at  0.5  sec  intervals.  The  peaks  should  and  do  occur  at  the  points 
±18.6°,  ±39.5°,  and  ±72.7°  which  correspond  (through  Equation  1)  with  the  0.5,  1, 
and  1.5  sec  delays.  It  is  also  seen  that  all  of  the  sampling  schemes  have  near¬ 
singularities  at  0°. 


Figure  16:  Condition  Number  of  Covariance  vs.  Bearing  Angle 
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4.4  Estimation  of  Bearing  Angle 

In  this  section  the  methods  of  estimation  of  bearing  angle  are  discussed.  The  two 
methods  used  in  the  estimation  of  bearing  are  the  minimum  mean  square  error 
(MMSE)  estimate  and  the  maximum  a  posteriori  estimate  (MAP).  The  data  used  in 
the  calculation  of  the  density  functions  used  in  the  calculations  to  follow  is  produced 
by  the  simulation  discussed  in  Section  4.2. 

The  first  method  of  estimation  that  is  considered  is  the  minimum  mean  square  error 
method.  The  density  function  for  bearing  is  calculated  for  source  bearings  ranging 
from  0 "  to  90".  For  each  bearing,  the  density  is  calculated  for  each  of  400  data  sets 
randomly  selected  from  a  large  simulated  data  set.  For  each  of  these  400  density 
functions  the  mean  is  calculated  using: 


•  \  (35) 

The  means  which  are  calculated  from  (1)  for  each  data  set  are  then  averaged  over  all 
data  sets  to  arrive  at  the  average  estimate  of  bearing. 


MMSE  -  -L  Ep, 
dsj-i 1 


(40) 


where  ds  is  the  number  of  data  sets. 
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The  second  method  employed  is  the  MAP  estimate.  In  this  method  the  same  400  data 
sets  are  used  to  calculate  the  densities.  In  this  method  the  estimate  is  made  by 
choosing  the  angle  which  corresponds  to  the  largest  value  of  the  density.  A 
maximum  is  found  for  each  of  the  400  density  functions  and  the  corresponding 
ffiimatw  of  angle  are  averaged  to  arrive  at  the  average  MAP  estimate. 


MAP  -  -L  ZMax,  <41> 

dsj-i  1 

where  Maxj  is  the  maximum  density  for  the  j'th  data  set.  Again,  ds  is  the  number  of 
data  sets. 

Figures  17,  18,  and  19  plot  the  MMSE  and  MAP  estimates  as  a  function  of  actual 
bearing  angle  for  various  sampling  schemes.  Figure  17  gives  the  estimates  for  a 
sampling  scheme  of  10  samples/sensor  at  1  sec.  intervals.  Figure  18  gives  the 
estimates  for  a  sampling  scheme  of  20  samples/sensor  at  1  sec.  intervals.  Figure  19 
gives  the  estimates  for  a  sampling  scheme  of  20  samples/sensor  at  0.5  sec.  intervals. 
The  sensor  separation  for  all  examples  is  X/4.  The  SNR  is  5.  The  confidence 
intervals  shown  in  Figures  17  through  22  have  90%  coverage  (see  Appendix  b). 
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Bearing  Angle  (deg) 

Figure  17:  Estimate  vs.  Actual  Bearing(10  Samples/  Sensor  -  2  sec  Interval) 
SNR-5. 


Bearing  Angle  (deg) 


Figure  18:  Estimate  vs.  Actual  Bearing  (20  Samples/Sensor  1  sec  Interval) 
SNR-5. 


57 


Bearing  Angle  (deg) 

Figure  19:  Estimate  vs.  Actual  Bearing  (20  Samples/Sensor  0.5  sec  Interval) 

SNR =5. 

It  is  seen  that  as  the  number  of  samples  increases,  the  estimates  of  bearing  more 
closely  reflect  the  actual  bearing.  This  is  also  true  for  the  case  of  increasing  the 
resolution  of  the  sampling  as  in  Figure  19.  In  all  cases  the  estimates  cross  each  other 
in  the  50*  -  60°  range.  This  indicates  that  the  density  functions  become  skewed  as  the 
actual  bearing  approaches  ±90°.  The  meaning  of  ‘skewed*  is  that  the  mean  value  of 
the  density  function  and  the  angle  which  corresponds  to  maximum  value  of  the  density 
function  are  not  the  same.  The  fact  that  the  MMSE  estimates  become  lower  than  the 
MAPs  indicates  that  the  mean  values  of  bearing  are  on  average  lower  than  the  values 
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of  bearing  for  which  the  density  is  a  maximum.  This  then  is  the  statistical  support  of 
the  observation  that  the  densities  become  skewed  as  bearing  angle  approaches  ±90° 
made  in  Section  3.5.  To  better  understand  these  trends  is  necessary  to  look  at  more 
statistics  derived  from  the  density  functions. 

4.5  Statistics  on  the  Estimates 

In  this  section,  the  statistics  of  the  estimates  are  presented.  In  the  last  section,  the 
average  estimates  of  bearing  are  calculated  using  a  large  number  of  data  sets  to 
calculate  the  conditional  probability  density  function  for  each  data  set.  Using  the 
same  techniques  as  in  the  previous  section,  the  variance  of  the  estimates  is  calculated. 
The  variances  of  the  MAP  and  MMSE  estimates  are  plotted  as  a  function  of  bearing 
together  with  the  Cramer-Rao  bound  on  the  variances.  The  statistics  are  different 
because  the  assumptions  about  the  bearing  made  in  the  thesis  are  fundamentally 
different  from  those  made  in  the  derivation  of  the  Cramer-Rao  bound  (C-R  Bound). 
The  C-R  bound  assumes  that  sufficient  prior  information  in  known  about  the  bearing 
so  that  the  covariance  of  the  data  is  known.  Under  this  assumption,  the  density  py,* 
is  Gaussian,  and  both  the  MAP  and  MMSE  estimates  are  identical.  When  the  bearing 
can  be  resolved  with  sufficient  accuracy  that  the  covariance  matrix  can  be  assumed 
constant  over  the  region  for  which  p^(Y  is  non-trivial,  the  C-R  bound  accurately 
reflects  the  performance  of  the  Bayes  optimal  estimators.  When  the  variations  of  the 
covariance  matrix  cannot  be  ignored  (as  happens  at  low  SNR)  the  C-R  bound  does 
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not  reflect  the  likely  performance  of  the  optimal  estimators.  It  is  important  to  note 
that  under  the  assumptions  used  to  derive  the  optimal  MAP  and  MMSE  estimators, 
the  former  estimate  is  not  unbiased,  and  neither  is  consistent,  therefore  both  may  have 
variances  lower  than  the  C-R  bound. 

The  Cramer-Rao  bound  is  the  lower  bound  on  the  estimate  variance  and  is  given  by 
the  relation 


Kor[*(fl)-i|r] 


d*2 


(42) 


where  \p  is  the  estimate  of  the  real  non-random  parameter  \p  and  Py^(Yj^)  is  the 
conditional  probability  density  function  of  the  data  r  conditioned  on  the  real  variable 
*• 

The  calculation  of  the  Cramer-Rao  bound  is  as  follows. 


The  logarithm  of  the  density  function  (4)  is  taken 

—  -  i 

The  second  derivative  with  respect  to  ^  is  taken 


(43) 


60 


-a  i!  i  <9  | 

r  to((2n  2  l*rl2)--^--(4l'r*r'‘n 
at2  3*2  2 


The  second  derivative  is  found  by  fitting  a  second  order  polynomial  to  (43)  near  and 
including  the  point  at  which  the  second  derivative  is  to  be  calculated.  This  is 
accomplished  by  first  finding  the  interpolating  polynomials  given  by 


//*)  =  D  *=1,2,3  (45 

A  V*y 

The  li’s  are  second  order  polynomials.  The  interpolation  requires  three  points:  the 
point  at  which  the  derivative  is  to  be  calculated  and  two  points  close  to  either  side. 
The  interpolation  function  is  given  by 


Mf)  -  2  xklkW 

k-Q 

where  the  Xk’s  are  coefficients  to  the  polynomials  given  by  the  functions 
ln((2x)<“/2)  1 RY  |  *)  and  (ViYTRy'Y)  evaluated  at  the  three  points  *2,  and 


The  second  derivative  for  a  second  order  polynomial  is  given  by  the  following 


*(♦) 


2*,* 


2*2* 


2*3* 


(*i-*2)(*»-*3)  (*s'*i)(*3-*2) 


(47) 
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The  Cramer-Rao  bound  is  found  by  taking  the  inverse  of  the  expectation  of  the 
relation  (43). 


°c* 1 


1 


£T 


-o  2.  1 


32 


(48) 


The  expectation  of  the  first  term  of  (43)  is  a  scalar.  Finding  the  expectation  of  the 
second  term  is  more  involved.  The  second  derivative  operator  is  moved  to  the 
inverse  of  the  covariance  giving 


R"  can  be  represented  by  a  singular  value  decomposition  by 

i 

R"-UH  Ut~QtQ 


(49) 


(50) 


A  new  variable  is  defined  as  Z=QY.  Using  this  new  variable  the  relation  for  the 
expected  value  of  the  second  term  on  the  right  hand  side  of  (43)  can  be  rewritten 
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E[YtR"  Y\  =  E[YtQ  tQY\=E[ZtZ]  =  Trace[E[ZZ  *))  (51) 

by  replacing  Z  with  the  proper  values  (51)  can  be  reduced  to 


%YtR"Y\  =  Trace[E[QYYTQT\]  *  TracelQ ElYY^Q7)  (52) 

The  term  on  die  far  right  of  the  (52)  contains  the  term  E[Y  Y7).  This  is  simply  the 
covariance  of  the  data.  Therefore,  the  entire  Cramer-Rao  bound  can  be  written  as 


1 


-  Jp  ln«2i.) W>  - 1 


(53) 


where  the  derivatives  are  calculated  as  mentioned  above  for  the  calculations  of  the 
first  term  on  the  right  hand  side  and  to  find  the  value  of  Q.  It  is  important  to  note 
that  (53)  is  for  unbiased  estimators.  In  our  case  the  variable  is  bearing  angle.  Note 
that  the  density  function  appearing  in  (4)  is  not  the  density  calculated  in  the  previous 
sections.  The  C-R  bound  is  interpreted  as  the  lower  bound  on  the  variance  of  the 
estimate  based  on  the  probability  density  of  the  data  conditioned  on  the  bearing  angle. 
It  is  important  to  keep  in  mind  that  the  density  function  calculated  in  all  of  the  MMSE 
and  MAP  estimates  is  the  density  function  of  the  bearing  angle  conditioned  on  the 
actual  data.  One  then  sees  that  the  C-R  bound  on  the  variance  is  arrived  at  by 
assuming  that  the  source  is  in  a  given  direction  and  then  based  on  the  data  gives  a 
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lower  bound  on  the  estimate  variance.  The  MAP  and  MMSE  estimates  are  made 
based  on  the  assumption  that  the  data  is  reflective  of  a  certain  bearing  angle.  In  other 
words  the  C-R  bound  assumes  the  source  is  in  a  given  direction  then  sees  how  close 
the  data  reflects  that  bearing.  The  method  used  to  calculate  the  variance  of  the 
estimates  in  the  MAP  and  MMSE  estimates  make  no  assumption  of  bearing  angle 
except  that  the  prior  density  on  bearing  is  uniformly  distributed.  The  estimates  are 
made  on  the  directly  calculated  probability  density  of  bearing  angle  based  on  the 
actual  data.  The  question  is  then  what  is  the  relationship  between  the  C-R  bound  and 
the  variance  of  the  estimates  using  our  calculated  density  functions.  This  question 
needs  to  be  asked  so  that  the  comparison  of  a  relatively  well  known  lower  bound,  the 
C-R  bound,  and  the  variances  of  the  estimates  made  using  the  directly  calculated 
density  functions  can  be  made.  The  answer  to  this  is  that  the  C-R  bound  shows  what 
land  of  results  are  obtained  using  an  efficient  estimator  such  as  a  maximum 
likelihood  estimator  using  the  known  covariance  of  the  data.  Figures  20,  21,  and  22 
show  the  MAP  and  MMSE  estimate  variances  along  with  the  C-R  bound. 
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Figure  20b:  Estimate  Variance  vs.  Bearing  (10  samples/Sensor  2  sec  Interval) 
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Figure  21a:  Estimate  Variance  vs.  Bearing  (20  Samples/Sensor  1  sec  Interval) 
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Figure  21b:  Estimate  Variance  vs.  Bearing  (20  Samples/Sensor  1  sec  Interval) 
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Figure  22a:  Estimate  Variance  vs.  Bearing  (20  Samples/Sensor  O.S  sec  Interval) 
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Figure  22b:  Estimate  Variance  vs.  Bearing  (20  Samples/Sensor  0.5  sec  Interval) 


67 


The  first  trend  that  is  seen  is  that  as  the  number  of  samples  increases  or  if  the 
sampling  resolution  increases,  the  C-R  bounds  get  lower.  It  is  also  seen  that  in  the 
10  sample  case  the  highest  variance  (worst  estimate)  is  at  a  bearing  angle  of 
approximately  65°.  This  is  an  unexpected  result.  One  would  expect  that  as  the  actual 
bearing  approached  ±90°  the  variance  of  the  estimate  would  increase  steadily.  To 
better  understand  why  this  occurs  one  has  to  look  at  the  averaged  density  functions. 

An  average  density  function  is  the  average  of  the  density  functions  which  have  been 
calculated  for  a  number  of  data  sets.  Figure  23  is  a  plot  of  the  averaged  density 
functions  for  400  data  sets.  One  sees  that  as  the  curves  approach  the  -65°  bearing 
angle,  the  densities  are  the  flattest  and  the  widest.  This  makes  the  estimation  of 
bearing  difficult  to  make  hence  the  largest  variance  on  the  estimate.  As  the  curves 
move  past  the  ±65°  bearing,  the  curves  become  sharper  with  higher  peaks.  This 
occurs  due  to  the  fact  that  the  density  functions  are  windowed  (in  the  range  ±90°)  by 
the  uniform  distribution  in  Equation  (18).  Since  the  curves  have  higher  peaks  and  are 
narrower,  the  variance  on  the  estimates  will  be  lower  in  this  range  and  continue  to  get 
lower  as  the  actual  bearing  of  the  source  approaches  ±90°. 
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Figure  23:  Averaged  Density  Functions  (10  Samples/Sensor  2  Sec  Interval) 

It  is  also  seen  that  low  points  or  "dips"  occur  in  the  C-R  curves  for  certain  sampling 
schemes  and  certain  bearings.  These  values  are  the  same  as  those  discussed  in  section 
4.3  for  which  the  condition  number  peaks.  It  is  at  the  values  which  the  covariance 
matrix  becomes  near  singular  that  the  variances  have  smaller  values.  By  comparing 
Figures  20,  21,  and  22  it  is  seen  that  the  near  singularities  corresponding  the 
sampling  interval  have  a  definite  effect  on  the  variances  of  the  estimates.  For 
example,  the  near  singularity  which  occurs  at  70*  in  Figure  22  for  a  sampling  scheme 
of  20  samples  /second  at  0.5  second  intervals,  brings  the  peak  in  the  variance  curve  at 
65°  down.  The  trend,  therefore,  in  this  particular  sampling  scheme  more  closely 
resembles  the  expected  trend  that  the  variance  of  the  estimates  of  bearing  goes  up  as 


bearing  angle  approaches  ±90*. 
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Looking  at  the  sampling  schemes  and  how  they  affect  the  variances  on  the  estimates  is 
also  interesting.  As  seen  in  Figure  20  the  MMSE  and  MAP  variances  are  smaller 
than  the  C-R  bound  for  bearings  greater  than  approximately  40°  and  for  the  MAP 
estimate  for  bearings  less  than  70°.  This  is  true  only  for  the  case  of  the  nyquist 
sampling.  The  MMSE  is  the  lowest  with  the  trend  being  that  there  is  a  peak 
occurring  around  30*  and  then  the  variances  drop  off  as  the  bearing  angle  approaches 
90*.  The  second  sampling  scheme,  where  20  samples  per  sensor  are  taken  at  1  sec. 
intervals  is  seen  in  Figure  21.  The  immediate  observation  is  that  the  variance  curves 
have  low  points  in  them  at  approximately  40°.  This  is  the  point  at  which  the 
condition  number  of  the  covariance  is  large  as  seen  in  Figure  16.  Since  the  variance 
here  is  the  lowest,  the  best  estimates  of  bearing  can  be  made  here.  It  is  also  noticed 
that  for  this  particular  sampling  scheme,  the  variances  are  all  closer  together.  The 
MMSE  is  still  generally  lower  than  the  MAP  and  C-R  bound,  but  the  C-R  bound  is 
now  clearly  lower  than  the  MAP  estimate  variance  curve.  Notice  that  a  peak  still 
exists  in  the  C-R  curve  at  approximately  63°.  The  reasons  for  this  are  that  same  as 
those  discussed  above  for  the  10  sample  per  sensoi  case.  Finally,  looking  at  the  20 
samples  per  sensor  case  at  0.3  second  intervals,  the  correlation  between  the  condition 
number  of  the  covariance  matrix  and  the  variance  plots  is  really  apparent.  All  of  the 
estimate  variance  curves  in  Figure  22  have  low  points  at  approximately  20*,  40*  and 
70"  bearings.  Looking  back  we  see  in  Figure  16  that  the  condition  number  plot  for 
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this  sampling  scheme  has  peaks  at  precisely  these  points.  Therefore  better  estimates 
can  be  made  for  a  sampling  schemes  with  more  samples  per  sensor  with  shorter 
intervals  between  samples.  This  is  true  due  to  the  fact  that  not  only  are  the  curves 
generally  lower,  but  singularities  occur  in  the  covariance  matrix  for  more  bearing 
angles  which  in  turn  causes  more  low  points  in  the  estimate  variance  curves. 


4.6  Summary  and  Conclusions 

In  this  chapter  a  new  simulation  was  created  with  a  more  complicated  noise  model. 
The  effect  of  this  new  data  model  on  the  covariance  matrix  of  the  data  is  discussed 
and  the  net  effect  is  found  to  be  that  block  diagonal  components  which  are  a  multiple 
of  the  signal  blocks  are  added  to  the  covariance.  This  is  true  since  the  noises  at  each 
sensor  are  correlated  with  themselves,  but  are  assumed  to  have  zero  cross-correlation 
between  sensors. 

Using  the  new  simulation,  new  data  is  created  which  is  used  to  calculate  the 
probability  density  functions  conditioned  on  this  data.  The  density  functions  are 
calculated  for  a  number  of  examples  using  various  sampling  schemes  on  the  data. 

The  density  functions  are  then  used  to  make  estimates  on  the  actual  bearing  angle. 
Two  types  of  estimates  are  made.  The  first  is  the  minimum  mean  square  error 
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estimate  (MMSE).  The  second  is  the  maximum  a  posteriori  estimate  (MAP).  It  is 
found  that  as  the  number  of  samples  is  increased,  or  if  the  sampling  interval  is 
shortened,  the  estimates  more  closely  reflect  the  actual  bearing. 

A  measure  of  how  accurate  the  estimates  are  is  then  discusses  in  the  form  of  the 
variances  on  the  estimates.  A  known  measure  of  performance  is  then  introduced:  the 
Cramer-Rao  bound  (C-R).  The  C-R  bound  gives  the  lower  bound  on  the  variance  of 
estimates  given  by  an  unbiased  estimator.  The  variances  on  the  MAP  and  MMSE 
estimates  are  lower  than  the  C-R  bound  for  Nyquist  time  and  spatial  sampling,  but  the 
C-R  bound  is  lower  for  any  form  of  over  sampling. 

Another  trend  is  found  in  relation  to  the  variances  on  the  estimates  and  the  condition 
number  of  the  covariance  matrix.  The  sampling  schemes  which  have  over-sampling 
contain  low  points  in  the  estimate  variance  curves  which  correspond  to  bearings 
where  the  covariance  has  high  condition  numbers  (covariance  matrix  approaching 
singularity).  For  Nyquist  sampling,  there  are  no  singularities  except  at  0°  bearing. 

As  the  sampling  interval  or  number  of  samples  is  significantly  increased,  the  number 
of  bearings  at  which  the  covariance  is  singular  '  1  increase  and  the  variances  on  the 
estimates  will  be  lower  which  indicates  that  better  estimates  on  the  actual  bearing  can 


be  made. 
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Chapter  5:  Summary  and  Conclusions 


5.1  Introduction 

In  this  chapter  the  results  obtained  in  the  thesis  are  discussed.  The  general  outline  of 
the  chapter  will  follow  the  outline  of  the  thesis  giving  the  results  obtained  thought. 
The  results  section  will  give  a  brief  description  of  the  goals  for  major  topics  in  the 
thesis.  The  opening  discussion  of  a  general  topic  is  followed  by  the  results  which  are 
presented  in  the  thesis  and  how  well  the  goals  of  the  section  are  fulfilled.  The 
chapter  ends  with  a  section  which  is  dedicated  to  recommendations  for  further 
research. 


5.2  Conclusions  and  Results 

In  Chapter  2  the  derivation  of  the  probability  density  function  conditioned  on  the 
sensor  inputs  is  presented.  It  is  shown  that  the  conditional  density  function  is 
calculated  by  first  calculating  the  probability  density  function  of  the  data  conditioned 
on  the  bearing  angle.  This  density  function  is  used  in  Bayes’  formula  to  calculate  the 
probability  density  function  of  the  bearing  angle  conditioned  on  the  sampled  data. 
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Following  the  derivation  of  the  probability  density  function,  an  analysis  of  the 
covariance  matrix  of  the  data  is  presented.  This  analysis  begins  with  the  placement  of 
the  data  into  an  ordered  vector  which  is  multiplied  by  its  transpose.  The  expected 
value  of  this  product  is  then  taken  to  give  the  covariance.  Following  this  discussion, 
the  code  which  is  used  to  generate  the  covariance  is  discussed  is  detail. 

In  Chapter  3,  a  simulation  used  to  generate  some  of  the  data  used  in  the  thesis  is 
presented.  The  elements  of  the  simulation  ore  discussed,  and  criteria  are  set  to  insure 
that  the  outputs  reflect  the  sensor  model.  The  outputs  of  the  simulation  are  analyzed 
in  order  to  give  meaning  to  definitions  of  signal-to-noise  ratio  and  the  ’whiteness’  of 
the  noise  processes. 

The  next  part  of  the  thesis,  also  in  Chapter  3,  are  four  sections  each  devoted  to  a 
different  example.  Each  section  contains  an  example  in  which  various  parameters  are 
changed.  The  parameters  changed  are  the  bearing  angle,  the  sensor  separation 
distance,  the  number  of  samples  taken  at  each  sensor,  and  the  signal-to-noise  ratio. 

In  the  first  example  it  is  found  that  as  bearing  angle  approaches  ±90°,  the  density 
functions  become  broader  and  have  lower  peaks.  This  occurs  up  to  a  point  at  which 
the  densities  begin  to  narrow  and  have  wider  peaks.  This  is  due  to  the  fact  that  the 
densities  are  cut  off  at  ±90°  .  The  second  example  shows  that  as  the  sensor 
separation  distance  is  increased,  the  number  of  possible  bearings  also  increases.  This 
is  due  to  spatial  aliasing.  The  separation  distance  is  larger  that  the  half-wavelength  of 
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the  center-frequency  of  the  signal  causing  a  sort  of  aliasing.  The  optimum  distance 
where  no  aliasing  occurs  is  the  Nyquist  separation  distance  where  the  separation 
distance  is  half  the  wavelength  of  the  signal  center  frequency.  In  the  third  example  it 
is  shown  that  as  the  number  of  samples  increases,  the  density  curves  become  more 
narrow  with  higher  peaks.  The  result  is  that  with  more  samples,  a  better  estimate  of 
bearing  angle  is  achieved.  The  fourth  example,  where  the  signal-to-noise  ratio  is 
changed,  shows  that  as  the  signal-to-noise  ratio  is  increased,  the  curves  become  more 
narrow  with  higher  peaks.  This  also  shows  that  as  the  SNR  is  increased,  better 
estimates  of  bearing  can  be  made. 

In  Chapter  4,  a  modification  of  the  simulation  discussed  in  Chapter  3  is  presented. 
This  new  simulation  has  broad-band  noise  exterior  to  the  narrow  band-pass  filters  in 
addition  to  the  broadband  signal.  The  noises  at  each  of  the  sensors  in  uncorrelated. 
The  net  effect  of  adding  this  noise  is  apparent  in  the  discussion  of  the  covariance. 
The  autocorrelation  of  the  noise  sequence  is  the  same  in  shape  as  the  autocorrelation 
of  the  signal,  but  its  magnitude  is  scaled  by  the  signal-to-noise  ratio.  The  effect  this 
has  on  the  covariance  is  that  block  diagonal  components  are  added  to  the  covariance 
which  are  scaled  versions  of  the  blocks  given  by  the  covariance  of  the  signal  alone. 
Only  the  diagonal  blocks  exist  since  the  noises  are  uncorrelated  between  sensors. 
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It  is  found  that  die  covariance  becomes  almost  singular  for  certain  bearings.  This 
occurs  for  bearing  angles  for  which  the  signal  propagation  delay  between  sensors  is 
equal  to  or  an  integer  multiple  of  the  sampling  interval. 

Chapter  4  moves  on  into  the  area  of  bearing  estimation  using  various  estimation 
techniques.  The  two  estimates  that  are  discussed  are  the  minimum  mean  square  error 
method  (MMSE)  and  the  maximum  a  posteriori  estimation  (MAP)  method.  These 
estimates  are  made  for  an  array  of  two  sensors  with  quarter  wavelength  spacing 
between  sensors  for  various  sampling  schemes.  It  is  found  that  for  a  large  number  of 
data  sets  (400)  the  MAP  and  MMSE  estimates  are  close  to  one  another  and  that  as  the 
number  of  samples  increases,  or  as  the  frequency  of  samples  increases,  the  estimates 
become  more  reflective  of  the  actual  direction  of  arrival.  It  is  also  seen  that  the 
curves  for  the  MAP  and  MMSE  estimates  intersect  at  a  point  around  ±60°.  This 
indicates  that  the  density  functions  become  skewed  and  are  no  longer  close  to 
Gaussian. 

In  order  to  better  understand  the  trends  seen  in  the  estimate  curves,  other  statistics 
such  as  the  variances  on  the  estimates  are  investigated.  These  variances  are  shown  to 
have  various  trends  according  to  the  sampling  scheme  used  to  calculate  the  density 
functions.  The  sampling  scheme  referred  to  Nyquist  sampling  has  differing  trends  for 
various  estimates.  The  trend  for  MAP  estimation  is  that  the  variance  curve  as  a 
function  of  bearing  increases  up  until  about  60°.  The  variances  then  get  smaller  until 
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±90°.  This  occurs  since  the  densities  are  limited  by  a  uniform  distribution  multiplier 
which  effectively  limits  the  densities  between  the  ranges  of  ±x/2.  The  MMSE 
curves  have  a  peak  in  their  variance  curve  at  approximately  ±30°.  The  variances 
then  drop  off  as  the  bearing  approaches  ±90°.  As  the  number  of  samples  is  increased 
and  the  sampling  interval  is  decreased,  low  points  or  'dips’  occur  in  the  variance 
curves  for  both  the  MMSE  and  MAP  estimate  variance  curves.  These  'dips'  occur  at 
exactly  the  points  mentioned  above  as  the  points  where  the  covariance  matrices 
become  singular. 

The  next  measure  of  the  accuracy  of  the  estimates  is  the  Cramer-Rao  lower  bound  on 
the  variance  of  estimates.  The  Cramer-Rao  bound  (C-R)  is  introduces  as  the 
variance  of  an  estimator  which  is  optimum  and  unbiased.  The  C-R  bound  is  used  to 
compare  the  estimates  made  using  the  density  functions  calculated  in  the  thesis  to 
estimates  which  would  be  made  by  an  optimum  estimator  using  the  same  covariance 
of  the  data.  The  trends  which  appear  in  the  plots  of  the  C-R  bound  vs.  bearing  angle 
are  similar  to  those  which  appear  in  the  MMSE  and  MAP  estimate  variances  which 
appear  above.  In  the  case  of  Nyquist  sampling,  the  trend  is  that  the  C-R  bound  as  a 
function  of  bearing  steadily  increases  until  approximately  ±65°.  The  curve  then 
drops  off  until  90°.  As  the  number  of  samples  or  frequency  of  sampling  is  increased, 
again,  'dips’  occur  in  the  C-R  lower  bound  curves.  These  'dips’  occur  precisely 
those  points  mentioned  above,  where  the  covariance  of  the  data  becomes  almost 
singular. 
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5.3  Recommendations  for  Further  Research 

In  this  section  a  few  ideas  for  further  research  are  presented.  The  problem  of 
direction  of  arrival  of  a  signal  relative  to  an  array  of  sensors  has  had  many  proposed 
solutions.  Many  algorithms  have  been  presents  for  the  solution  of  the  problem.  The 
object  of  this  section  is  not  to  suggest  new  algorithms,  but  how  to  use  the  results  of 
the  thesis  in  order  to  test  the  performance  of  algorithms  which  already  exist.  The 
C-R  bound  is  often  used  as  a  performance  test  by  comparing  the  estimate  variance  to 
the  C-R  lower  bound.  The  estimates  on  the  bearing  angle  made  in  this  thesis  are 
done  using  a  different  density  function  from  that  used  in  the  C-R  bound  tests.  The 
density  functions  calculated  in  this  work  are  conditioned  on  the  actual  data  received  at 
by  the  array.  The  C-R  densities  used  in  calculating  the  bound  are  conditioned  on  an 
assumed  value  of  bearing  angle.  The  difference  is  subtle,  but  the  result  is  that  the 
estimate  variances  presented  is  Chapter  4  are  an  alternative  test  of  performance. 


The  first  recommendation  is  that  a  code  be  written  which  is  efficient.  The  code 
presented  in  this  thesis  which  computes  the  conditional  probability  density  function 
has  a  very  long  run  time.  For  each  point  in  the  curves  which  are  presented  in 
Chapter  4,  a  large  quantity  of  CPU  time  is  required.  An  efficient  code  in  FORTRAN 
or  ADA  could  significantly  reduce  the  run  time.  If  the  run  time  were  significantly 
reduced,  runs  could  be  d'Mie  to  calculate  density  functions  conditioned  on  larger  data 
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sets.  In  this  thesis  the  largest  data  set  consists  of  twenty  points.  It  would  by 
interesting  to  see  what  effect  a  large  number  of  near-singularities  occurring  in  the 
covariance  matrix  has  on  estimates  of  bearing  and  the  corresponding  variances  on 
those  estimates. 

The  next  interesting  area  of  research  would  be  to  use  the  algorithm  developed  to 
calculate  the  conditional  density  function  of  the  bearing  angle  for  real  data  rather  than 
simulated  data.  Some  modifications  to  the  code  may  me  necessary  to  include  possible 
noise  correlation  effects  between  sensors.  The  algorithm  should  also  include  more 
than  two  sensors.  The  work  done  above  used  only  two  sensors  due  to  limited 
computing  resources.  It  would  be  interesting  to  compute  the  density  functions  for  real 
data  for  numerous  sensors  and  compare  the  estimates  of  bearing  angle  made  using  the 
calculated  densities  vs.  the  bearing  estimated  using  spectral  or  beamforming 
techniques  discussed  in  Chapter  I. 


Another  interesting  area  of  research  would  be  to  see  what  effect  random  sampling 
intervals  would  have  on  the  density  functions  and  the  estimate  variances.  The  effects 
on  the  covariance  singularities  and  the  related  effects  on  the  estimates  would  be 
interesting  to  see. 
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Appendix  A 


The  following  is  a  listing  of  the  code  which  generates  the  subblocks  of  the  covariance 
matrix  discussed  in  Chapter  2. 


•♦The  subroutine  ’User’  decides  which  subblock  is  to  be  filled  then  calles  the  necessary 
subroutines  which  fill  the  appropriate  subblock. 

subroutine  user(a,ma,g,s,t) 

real*8  a(1025, 100),dpn,tau0,s,t,rs0(1025,  l),tau(1025, 1) 

integer  n,pn,sens,sensl,12,13,ma,g 
n=a(l,3) 
pn«a(2,3) 
dpn=a(3,3) 
tau0=a(4,3) 
n0=a(5,3) 
sens*a(6,3) 
sensl  =a(7,3) 
do  2  dd- 1,1025,1 
tau(dd,l)=a(dd,2) 
rsO(dd,l)==a(dd,l) 
a(l,l)«rs0(513,l) 
a(2,2)=rs0(513,l) 
a(3,3)=rs0(513,l) 

2  continue 

if  (sens.eq.sensl)  then 

call  diagonal^, pn,dpn,d,a,rsO,tau,nO) 

end  if 

if  (sens.lt.sensl)  then 
do  40 12-l,pn,l 
do  30 13  **  1  ,pn,  1 
if  (12.eq.13)  then 

call  wowl  (12,13,tau0,sens,sensl,a,dpn,rs0,tau,n0) 
end  if 

if  (12.lt.13)  then 

call  wow2  (12,13,tau0,sens,sensl,a,dpn,rs0,tau,n0) 
end  if 

if  (12.gt.l3)  then 

call  wow3  (12, 13, tauO, sens, sensl, a,dpn,rsO,tau,nO) 
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end  if 

50  continue 

40  continue 

end  if 

if  (sens.gt.sensl)  then 
do  60  12=l,pn,l 
do  70  13=l,pn,l 
if  (12.eq.I3)  then 

call  wow4  (12,13, tauO, sens, sensl, a, dpn,rsO,tau,nO) 
end  if 

if  (12.lt.13)  then 

call  wow5  (12,13,tau0,sens,sensl,a,dpn,rs0,tau,n0) 
end  if 

if  (12.gt.13)  then 

call  wow6  (12,13,tau0,sens,sensl,a,dpn,rs0,tau,n0) 
end  if 

70  continue 

60  continue 

end  if 
return 
end 


**  The  subroutine  ’Wowl*  fills  the  main  diagonals  of  the  subblocks  which  lie  below 
the  main  diagonal  subblocks. 

subroutine  wowl(12, 13, tauO, sens, sensl ,a,dpn,rs0,tau,n0) 
real*8  block(200, 200), tauO, a(1025, 100), dpn,rs0(1025, 1), 

#tau(1025,l) 

integer  n0,sens, sensl  JlJhJmJ, 12,13 
block(13 ,12) = (sens  1  -sens)  ♦tauO 
jh=n0 
jl*l 

jm=int((jh-jl)/2) 
do  While  ((jh-jl)  .gt.  1) 

if  (block(13,12).gt.tau(jm,l))  then 
jl=jm 

jm=int(jm+((jh-jl)/2)) 

else 

jh=jm 

jm  *int(jm-((jh-jl)/2)) 
end  if 
end  do 
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j=jl 

a(13,12)=((RsO(j  + 1 ,  l)-rsO(j,l))/(tau(j  + 1 ,  l)-tau(j ,  1)))* 
#(block(13,12)-tau(j,  l))+rsO(j,  1) 
return 
end 


**  The  subroutine  ’Wow2’  fills  the  area  below  the  main  diagonal  of  the  subblocks  which 
lie  below  the  main  diagonal  subblocks. 

subroutine  wow2(12,13,tau0,sens,sensl  ,a,dpn,rsO,tau,nO) 
real*8  a(1025, 100),block(200,200),tau0,dpn,rs0(1025, 1), 

#tau(1025,l) 

integer  sens,sensl,12,13,n0jhjljmj 
block(13 , 12) = ((sens  1  -sens)  *tauO) + ((13-12)  *dpn) 
jh*nO 
jl=l 

jm=int((jh-jl)/2) 
do  WhUe  ((jh-jl)  ,gt.  1) 

if  (blockG3,12)  .gt.  tau(jm,l))  then 
jl«jm 

jm  *int(jm +((jh-jl)/2)) 

else 

jh=jm 

jm=int(jm-((jh-jl)/2)) 
end  if 
end  do 

j=jl 

a(13,12) » ((RsO(j + 1 ,  l)-rsO(j ,  1))/ (tau(j + 1 ,  l)-tau(j  ,1)))* 

#(block(13,12)-tau(j ,  1)) +rsO(j ,  1) 
return 
end 
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**  The  subroutine  ’Wow3’  fills  the  area  above  the  main  diagonal  of  the  subblocks  which 
lie  below  the  main  diagonal  subblocks. 


subroutine  wow3(12, 13, tauO, sens, sens l,a,dpn,r$0,tau,n0) 
real*8  block(200,200),tau0,a(1025,100),dpn,rs0(1025,l), 
#tau(102S,  1) 

integer  sens,sensl,12,13,n0jh jm,j  jl 
block(13,12) = ((sens  1  -sens)*tau0)-((12-13)*dpn) 
jh=nO 
jl=l 

jm=int((jh-jl)/2) 
do  While  ((jh-jl)  .gt.  1) 

if  (block(13,12)  .gt.  tau(jm,l»  then 
jl=jm 

jm*int(jm+((jh-jl)/2)) 

else 


jh=jm 

jm=*int(jm-((jh-jl)/2)> 
end  if 
end  do 


j 

a(13,12) =((RsO(j + 1 ,  l)-rsO(j,  1))/ (tau(j + 1 ,  l)-tau(j ,  1)))* 
#(block03,12)-tau(j,  l))+rsO(j,  1) 

return 

end 


**  The  subroutine  ’Wow4’  fills  the  main  diagonals  of  the  subblocks  which  lie  above 
the  main  diagonal  subblocks. 


subroutine  wow4(12,13,tau0,sens,sensl  ,a,dpn,rsO,tau,nO) 
real*8  block(200,200),tau0,a(1025, 100),dpn,rs0(1025,l), 
#tau(1025,l) 

integer  sens,sensl,12,13,n0jhjljmj 
block(13 ,12) = ((sens-sens  1 )  *tauO) 
jh=nO 
jl«l 

jm=int((jh-jl)/2) 
do  While  ((jh-jl)  .gt.  1) 

if  (block(13,12)  .gt.  tau(jm,l))  then 
jl=jm 

j  m = int(j  m + ((j  h-ji)/2)) 
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else 

jh=jm 

jm=int(jm-((jh-jl)/2)) 
end  if 
end  do 

j=jl 

a(12,13) = ((RsO(j  + 1 ,  l)-rsO0 ,  l))/(tau(j  + 1 ,  l)-tauG ,  1)))* 
#(block(13 ,12)-tau(j ,  l))+rs0fl ,  1) 
return 
end 


**  The  subroutine  ’Wow5*  fills  the  area  below  the  main  diagonal  of  the  subblocks  which 
lie  below  the  main  diagonal  subblocks. 


subroutine  wow5(12,13,tau0,sens,sensl  ,a,dpn,rsO,tau,nO) 
real*8  block(200, 200),tau0,a( 1025, 100), dpn,rsO( 1025,1), 
#tau(1025,l) 

integer  sens,sensl,12,13,n0jh,jljmj 
block(13,12) = ((sens-sensl)*tau0)+((13-12)*dpn) 
jh=n0 
jl=l 

jm=int(Gh-jl)/2) 
do  While  ((jh-jl)  .gt.  1) 

if  (block(13,12)  .gt.  tauGm,l))  then 
jl=jm 

j  m = intG  m + ((j  h-j  l)/2)) 

else 

jh=jm 

jm =int0'm-(0‘h-jl)/2)) 
end  if 
end  do 

j-jl 

a(12,13)=((RsO(j + 1 ,  l)-rsO(j ,  l))/(tauG + 1 ,  l)-tau(j ,  1)))* 
#(block(13,12)-tauG,l))+rsO(|,l) 
return 
end 
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**  The  subroutine  ’Wow6’  fills  the  area  below  the  main  diagonal  of  the  subblocks  which 
lie  below  the  main  diagonal  subblocks. 

subroutine  wow6(12,13,tau0,sens,sensl  ,a,dpn,rsO,tau,nO) 
real*8  block(200,200),tau0,a(1025,100),dpn,rs0(1025,l), 

#tau(1025,l) 

integer  sens,sensl,12,13,n0jh jljmj 
block(13,12) = ((sens-sens  l)*tau0) + ((13-12)*dpn) 
jh=nO 
jl*l 

jm=int((jh-jl)/2) 
do  While  (Cjh-jl)  -gt.  1) 

if  (block(13,12)  .gt.  tau(jm,l»  then 
jl=jm 

jm =int(jm +((jh-jl)/2)) 

else 

jh*jm 

jm =int(jm-((jh-jl)/2)) 
end  if 
end  do 

j=jl 

a(12,13) =((RsO(j + 1 ,  l)-rsO(j,  l))/(tauO +1,1  )-tauQ,  1)»* 

#(block(13,12)-tau(j,  l))+rsO(j,  1) 
return 
end 


**  The  subroutine  ’Diagonal  nils  the  main  diagonal  subblocks. 


subroutine  diagonal(n  ,pn ,  dpn ,  d ,  a ,  rsO ,  tau ,  nO)  * 
real*8  delay(2, l),d(200,200),dpn,a(1025, 100), 
#rs0(1025,l),tau(1025,l) 
integer  n,pn,sn, shift, nOjh  jl,jmj,czap,rzap 
do  10  b_column=l,n,l 
do  20  b_row=l,n,l 
sn=0 

if  (b_row.eq.b_column)  then 
do  25  ll  =  l,pn,l 
Delay  (11 , 1) =01  *dpn)-dpn 
25  continue 

shift =b_column-l 

do  30  column =l,pn 


do  40  row=2,pn-sn 
czap = column + b_columi  l + shift 
rzap = b_row + row + sn + shift- 1 
d(rzap,czap) = delay  (row,  1) 
d(czap,rzap)  =delay(row,  1) 
jl=l 
jh«nO 

jm=int((jh-jl)/2) 
do  while  ((jh-jl).gt.l) 
if  (d(rzap,czap).gt.tau(jm,l))  then 
jl-jm 

jm=int(jm+((jh-jl)/2)) 

else 

jh=jm 

jm =int(jm-((jh-jl)/2)) 
end  if 
end  do 

a(rzap,czap) =((RsO(jl + 1 ,  l)-rsO(jl,  l))/(tau(jl+ 1 , 1)- 

tau(jl,l)))* 

#(d(rzap,czap)-tau(jl,  l))+rsO(jl,  1) 
a(czap,rzap) = ((RsO(jl +1,1  )-rsO(ji,  l))/(tau(jl +  1,1)- 
tau(jl,l)»* 

#(d(rzap,czap)-tau(jl,  1))+isOQ1,  1) 


40 

continue 

sn=sn+l 

30 

continue 

end  if 

20 

continue 

10 

continue 

return 

end 


Appendix  B:  Confidence  Intervals 


In  this  appendix,  the  method  used  to  calculate  the  confidence  intervals  used  in  chapter 
4  are  discussed.  The  confidence  intervals  are  used  in  Chapter  4  since  the  discussion 
covers  calculations  of  statistics  using  large  data  sets.  The  discussion  which  follows  gives 
a  brief  description  of  the  method  used  to  calculate  confidence  intervals  which  have  90% 
coverage  on  calculate  parameters. 

For  a  given  set  of  data  Xi,  X2,  X3,  X,,  with  a  finite  mean  n  and  sample  variance 

s*(n),  the  confidence  interval  based  on  the  t  distribution  for  large  n  is  given  by 

(» 

(  r  J  n 

The  parameter  a  is  determined  by  die  percent  confidence  interval  given  by  100(1 -a). 
The  parameter  is  the  upper  l-a/2  critical  point  for  a  standard  normal  distribution. 
See  table  T.  1  of  [1].  The  interpretation  of  the  confidence  interval  as  given  by  (1)  is  as 
given  by  [1] 


"If  one  constructs  a  very  large  number  of  100(1 -a)  percent  confidence  intervals  each 
based  on  n  observations,  where  n  is  sufficiendy  large,  the  proportion  of  these  confidence 
intervals  which  contain  (cover)  n  should  be  (1-a).  This  proportion  is  called  the  coverage 


for  the  confidence  interval." 


